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ABSTRACT 


vAo» N n W A and , surface hydrological parameterizations are implemented into the 
Institute for Space Studies (GISS) General Circulation Model 
(GCM). These parameterizations are: 1) runoff and ovapotranspiration 

functions that include the effects of subgrid scale spatial variability and use 
physically-based equations of hydrologic flux at the soil surface, and 2) a 
realistic soil moisture diffusion scheme for the movement of water in the soil 
column. 

A one-dimensional climate model with a complete hydrologic cycle is 
used to screen the basic sensitivities of the hydrological parameterizations 
before implementation into the full three-dimensional GCM. Results of the 
final simulation with the GISS GCM and the new landsurface hydrology 
indicate that the runoff rate, especially in the tropics, is significantly 
improved As a result, the remaining components of the heat and moisture 
balance show comparable improvements when compared to observations. 

The validation of model results is carried from the large global (ocean 
and landsurface) scale, to the zonal, continental, and finally the finer river 
basin scales. 


(Key words: Climate modeling, Global hydrology) 
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Chapter I 

Background and Introduction To Model Formulation 


A. General Circulation Models 

1. Brief Description 

The climatic system is exceedingly complex; the interaction of physical 
processes producing climate cover an enormous range in scale of time and 
space. The accurate scientific representation of a climate system in full detail 
is virtually impossible. There is no laboratory in which one may carry out 
controlled experiments on climate. Furthermore, the sampling density of most 
key climatic parameters turn out to be, for practical reasons, quite sparse over 
many climatic regions. Some parameters such as evaporation cannot even be 
measured directly. Nevertheless the scientific study of climate proceeds in 
spite of the inherent difficulties. 

In light of the controversial threat of global warming, what would most 
be desired is some means of understanding the effects of changes in the 
boundary conditions of the climate system. Today, the best analog we have 
to climate on a global scale, with which we can also perform experiments, is 
the General Circulation Model (GCM). The following is a very brief and 
abbreviated description of GCMs. For a fuller treatment, there are several 
texts available (e.g., Henderson-Sellers and McGuffie, 1987, and Washington 
and Parkinson, 1986). 

The first GCM was pioneered in 1956 by N. A. Phillips. This model 
had a simple 2-layer atmosphere over the Northern Hemisphere which 
incorporated quasi-geostrophy and hydrostatic equilibrium. A finite difference 
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scheme solved these equations over a 17 * 16 point gridded area covering 
60 million km 2 . The model succeeded in producing a jet stream and the 
3-celled structure observed in the earth’s atmosphere (the Hadley, Ferrel, and 
polar cells). Since then GCMs have increased in sophistication. Today for 
instance, numerical solution to the so-called primitive equations (see Table 1.1) 
has superseded the quasi-geostrophic assumption; the domain has become 
global instead of hemispheric; detailed radiative transfer schemes have been 
developed; more realistic boundaries (e.g., topography, albedos, emissivities, 
etc.) have been added, and atmospheric interactions producing precipitation 
have been refined. In the late 1960s, S. Manabe pioneered the effort to 
incorporate a hydrologic cycle. Before this, GCMs were "dry" and the latent 
heat sources and sinks were parameterized. 

There are several GCMs in existence which continue to be improved and 
refined as new techniques are implemented and as computer capabilities are 
expanded. Three of the major GCMs in the United States often referenced for 
their predictions regarding global warming are 1) the National Oceanic and 
Atmospheric Administrations’s (NOAA) Geophysical Fluid Dynamics Laboratory 
(GFDL) Model in Princeton, NJ, 2) the National Center for Atmospheric 
Research (NCAR) Community Climate Model (CCM) in Boulder, CO and 3) 
the National Aeronautic and Space Administration’s (NASA) Goddard Institute 
for Space Studies (GISS) Model in New York, NY. The GISS Model II (i.e., 
Version II, or 2nd generation, described in Hansen et al., 1983) has been the 
model used in this research. 

For GCMs to be of practical value, they need to be simplified (in terms 
of the resolution of discretization) in order to meet computational constraints. 
Thus, typical GCM grid spacing is of the order of hundreds of kilometers to a 
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Table 1.1 Fundamental ("Primitive") Equations of the atmosphere used in 
GCMs (after Hansen et al., 1983). The six unkn owns in the 
atmospheric state vector (/?, P, T, $) are numerically solved for, 
using Equations (T1)-(T6). 






Conservation of Momentum 

dV 

dt 

-4-4 -4 

= — 2fi x V - /f 1 V P + q + F 

(T1.T2.T3) 

Conservation of Mass 

d£ 

dt 

T > 
o 
II 

(T4) 

Conservation of Energy 

dU 

dt 

- - p 5f 1 + « 

(T5) 

Ideal Gas Law 

P - 

= pRT 

(T6) 


(Equation of State) 

-» 

where V = (llj + Uj + Wfc) velocity relative to rotating earth 

t time 

-4 

fi planet’s angular rotation vector 

p atmospheric density 

—4 — i 

g apparent gravity [== true gravity — Q x (fi * r)] 

r position relative to planet’s center 

-4 

F force per unit mass 

U internal energy per unit mass [= c^T] 

Q heating rate per unit mass 

R gas constant 

c v specific heat at constant volume 
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grid square side. The grid spacing of the GISS Model II for 8° * 10° grids 
(it can also be run using 4° * 5° grids) is shown in Fig. 1.1. At this 

resolution, running time on the IBM mainframe computer requires roughly one 
hour computing time per simulated month. 

GCM grid resolution allows solution of the primitive equations but 

typically is much too coarse for many of the other important processes. Thus, 
phenomena such as radiation, moist convection, condensation, cloud formation, 
boundary layer interaction, and landsurface hydrology must be "parameterized." 
This term describes the simplification of complex processes to accommodate 
computational constraints. The simplest of all parameterizations would be to 
leave out a process altogether. More commonly, a process is represented by 
simplified empirical relationships having dependence on a greatly reduced 

number of variables, often only one, as we shall see in the case of current 
hydrologic parameterizations. 

Because of the large amount of simplifications, and the coarse grid 

spacing, and the lack of understanding of important physical phenomena such 

as cloud formation, GCMs have been the subject of controversy in scientific 
circles, particularly in regard to global warming predictions. It has been 

pointed out (Stone and Risbey, 1990) that a key element critical to the 

improvement of GCMs continues to be the development of reliable subgrid 
scale parameterizations. The goal of this work is to show the vast 
improvements in the GISS landsurface hydrology to be realized by the 
implementation of the subgrid scale parameterization of Entekhabi and 

Eagleson (1989a). 
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Figure 1.1 


Grid discretization of the GISS GCM using 8° * 10° resolution. 


9 Improving the T.a.nH surface Hydrological Parame terization in the 


GTSS GCM: An Overview of this Research 

Precipitation in GCMs is produced, for each grid and at each time step, 
according to the physics and thermodynamics of the overlying air column. 
Similarly the potential evaporation is computed as the surface latent heat flux 
under conditions of unlimited moisture supply. The task of the landsurface 
hydrology routine is to partition the incident precipitation into runoff and 
infiltration and reduce the potential evaporation to the actual value. The 
partitioning and reducing rates depend on the state of the near-surface soil 
moisture; therefore the landsurface routine also tracks soil moisture balance. 

The current landsurface parameterization in the GISS GCM uses 
simplified empirical equations; they are chosen for their computational 
efficiency and adjusted to create reasonable gross estimates of global runoff 
while also maintaining reasonable surface temperatures over continental regions. 
» ... [the Model II runoff scheme] yielded 24 cm year -1 global runoff and 
summer continental temperatures about 2— 3°C cooler than observed. This 
intermediate runoff is used in Model II, primarily for its effect on summer 
temperatures", (Hansen et al., 1983). As will be discussed in detail in 
Chapter III, there are many shortcomings of this current hydrological 
parameterization, particularly in runoff and evaporation over tropical regions. 

As an improvement on GCM hydrological parameterizations in general, 
Entekhabi and Eagleson (1989a) developed new formulations for runoff and 
evaporation. The technique is based on reasonable spatial probability 
distributions of rainfall and soil moisture; the model updates these variables 
using the means (first moments) of the distributions. As is illustrated in 
Fig. 1.2, the probability density function (pdf) for precipitation intensity P in 
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Uniform vs. Distributed Precipitation 



Precipitation Intensity, P 


Figure 1.2 


Generalized probability distributions of distributed versus 
uniform precipitation intensity 



the GISS model (and all other GCMs used today) is a delta function at the 
model update (mean) precipitation value. That is, the precipitation is 
spatially uniform since all probability "mass" is concentrated at the one mean 
value. Eagleson and Entekhabi (1989a), on the other hand, incorporate a 
distributed pdf of precipitation which recognizes that rainfall, in nature, is very 
far removed from a simple uniform structure over the entire grid surface. 
This curve in Fig. 1.2 assumes that rainstorms consist of large areas of low 
intensity rainfall and relatively smaller areas of intense precipitation. Later on 
it will be shown that by assigning a specific mathematical function to this 
distribution and also one for soil moisture distribution, analytic expressions 
may be derived for runoff and evaporation. The resulting parameterization has 
the distinct critical advantage of accounting for spatial variability while being 
computationally efficient for practical use in GCMs. Preliminary testing of 
this new parameterization was done on a simplified one-dimensional climate 
model developed at MIT (Koster et al., 1988, and later modified by Entekhabi 
and Eagleson, 1989b) based on the GISS GCM. 

The goal of this work is to expand the current one-dimensional model to 
include such things as a soil heat diffusion mechanism (to accommodate 
seasonal simulations) and the inclusion of root extraction of soil moisture to 
simulate plant transpiration. Sensitivity will be performed by testing these 
various aspects in Chapter II. Finally, Chapter III will be devoted to the 
implementation of the new runoff and evaporation formulations in the GISS 
GCM. The new soil moisture diffusion scheme of Abramopoulos et al. (1988) 
will also be implemented with root extraction from lower layers. An 
evaluation of the results will be done using multiple observation data sets for 
the global hydrologic cycle and surface temperature. 
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3. One-Dimensional Modeling as an Efficient Screening Tool 

Long-duration simulation experiments with a full three-dimensional GCM 
for the purpose of testing landsurface hydrological parameterizations is an 
expensive proposition. It is therefore necessary to use a simpler and more 
efficient model to screen the many possible simulation experiments so that 
comprehensive testing in the full three-dimensional GCM is reserved for 
well-defined and focused experiments. In the case of the one-dimensional 
model at MIT, computer simulation time requirements are reduced from that 
of the GISS GCM by a factor of roughly 20. The reduced model must 
nonetheless contain the major interactions inherent in the full GCM or else the 
sensitivities in the reduced model and the full GCM will not agree with one 
another. 

The main compromise made to achieve such efficiency lies in the 
reduction to one dimension. Obviously with only one grid square there can be 
no network of nodal points over which one may solve the equations of 
atmospheric motion. Thus pressure levels and winds must be prescribed. The 
one-dimensional model of MIT, hereafter referred to as the Screening Model, 
dynamically links one ocean and one land grid in order to contain a complete 
hydrologic cycle. 

The dynamic exchange of heat and moisture between these two grid 
squares is parameterized. Having thus a simplified ocean source of heat and 
moisture, and recognizing that most of the physical processes such as radiation, 
moist convection, and condensation are one-dimensional (vertical in the 
atmospheric column) even in the GISS GCM, the Screening Model landsurface 
grid square retains sensitivities and interactions similar to the full 
three-dimensional GCM. As such it is an effective tool for landsurface 


21 



hydrological parameterization testing. The screening model will be described in 
greater detail in subsequent sections. 

B. Basic GCM Hydrological Parameterizations 

1. The Water Balance 

The purpose of hydrological parameterizations in GCMs is essentially to 
account for regional water balance. Each grid square variable is a point 
representative of a large area. Precipitation formed in the atmospheric column 
above a landsurface grid square is partitioned into two components, one 
infiltrating into the soil and one leaving the grid as surface runoff. Generally 
there is no lateral surface or subsurface flow between grid squares, and so 
water movement in the soil is restricted to the single vertical dimension. 
Surface runoff is eliminated instantaneously from the landsurface, whether 
adjacent to the ocean or not, under the assumption that it is efficiently 
delivered to the ocean via rivers. The final component of the terrestrial water 
balance is then the evapotranspiration from the surface of the soil column. 
Mathematically we have for each landsurface grid 

= P - e - Q ( ll > 

where 

W water stored in soil column (per unit area) 

P precipitation 

e evapotranspiration 

Q surface runoff 

In nature there is a wide variety of detailed processes interacting to 
produce the final water balance. For example, atmospheric heterogeneities 
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create highly irregular rainfall patterns spatially and temporally; soil texture 
heterogeneities in combination with topography result in wide variations in 
runoff and infiltration; myriad varieties, shapes, and sizes of vegetation 
influence evapotranspiration; fluctuations in the local heat budget cause 
variability in the ambient temperature, etc. These processes are simplified or 
neglected altogether in GCMs to accommodate computational constraints. 

Runoff and evapotranspiration typically are taken as fractions of 
precipitation and potential evaporation respectively; empirical relationships with 
dependence only on relative soil saturation (ratio of water volume to void 
volume; range: between 0 and 1) are used. Moisture diffusion in the soil is 

typically a simple mass balance relationship. Some more recent hydrological 
parameterizations are incorporating much greater detail and are 

physically-based, but very few offer any sophisticated means of accounting for 
spatial heterogeneity over the large extent of grid squares. The 
parameterization of Entekhabi and Eagleson (1989a) is an exception. 
Section I.B.2 will summarize two very simple parameterizations (those of GISS 
and GFDL); two of more detail will be described in Section I.B.3 (Sellers 
et al., 1986, and Dickinson et al. , 1986). Section I.C reviews the inclusion of 
subgrid scale spatial variability using the approach of Entekhabi and Eagleson 
1989a). 

2. Current GISS and GFDL Hydrology Schemes 

The storage component, W, in (1.1) is handled by a simplified diffusion 
between two layers of soil in the GISS GCM. There is also a surface snow 
storage during cold seasons. Our concern here, though, is mainly the 
formulations of runoff and evapotranspiration. 
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Runoff is taken to be a fraction R of the precipitation. The value 


Q = RP 0 < R < 1 ( 12 ) 

of R in the GISS Model is taken to be one-half the value of the relative 


surface soil saturation s 


R = 1/2 s 0 < s < 1 


(1.3) 


where s is defined as 


where W is again the per unit area stored soil water depth. Any infiltrating 
water which would exceed soil field capacity is also drained as runoff. 

Evaporation in a similar way is taken as a fraction 0 of potential 

evaporation e p . We have 

e = /?e p 0 < 0 < 1 ( L5 ) 

The value of 0 , like R, is parameterized as a simple function of relative 
surface soil moisture, 

0 = s 0 < s < 1 ( 16 ) 

In the GFDL Model (Manabe et al., 1969), the runoff ratio R is 
replaced by a "bucket" model in which runoff is only produced by an 
exceedance of soil field capacity. Therefore, for the GFDL model 


P - e W = W pc 

0 w < w FC 


where W pc is the soil field capacity. 


(1.7) 
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The GFDL evaporation function is of the same form as (1.5) but (3 is 
given by 


P = 



s. > s < 1 

k ~ 

0 < s < s. 

k 


( 1 . 8 ) 


Here s^ is a "critical value" of s usually taken to be 0.75. 

These very simple empirical relationships, very loosely based on 
hydrologic research, are used mainly for their computational ease while 
providing at least a gross similarity between model output and observed data. 
Two other parameterizations which incorporate a great deal more detail are 
the SiB (Simple Biosphere, Sellers et al., 1986), and BATS (Biosphere 
Atmosphere Transfer Scheme, Dickinson et al., 1986). 

3. SiB and BATS 

The parameterizations of SiB and BATS are improvements insofar as 
they allow for much greater detail in vegetation and the other surface 
processes along with more realistic soil diffusion. The basic feature of these 
two parameterizations is their representation of the soil-vegetation-atmosphere 
transfer of heat, moisture and momentum through resistance formulation. Here 
the flux of the property is proportional to the gradient in the property 
evaluated at two levels; the constant of proportionality is the 

inverse— resistance. The plethora of detail involving such things as leaf 
temperature and leaf angle, stomatal and root resistance, canopy heat capacity, 
interception capacity of leaves, and extinction coefficients of canopies, among 
others, is too broad to allow a concise summary of equations. The main 
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advantage is the representation of many processes which contribute to the 
hydrologic cycle as well as the radiation and heat balance. The difficulty of 
such an approach, however, is that most of the equations are still empirical 
and furthermore require numerous parameters which must be estimated from 
sparse data or approximated in some way where data is nonexistent. Also, 

there are stronger computational demands. 

The element missing from both the detailed SiB and BATS schemes as 
well as from the GFDL and GISS schemes, is the accounting for spatial 
variability across the large GCM grid area. 

C. Hydrological Parameterization with Spatial Variability 

The following is a condensed version of the derivation of the runoff 
coefficient R and evapotranspiration efficiency (3 after Entekhabi and Eagleson 
(1989a). 

1. Runoff Coefficient 

As with the GISS model, runoff Q is also taken as a fraction R of 
precipitation (see Eq. 1.2). It is through the formulation of R that the 
dependence on spatial variability is included. After Warrilow et al. (1986) 
precipitation intensity is assumed to be exponentially distributed over a wetted 
fraction k of the grid square, 



kP 
E [P] 


P > 0 


(1.9) 


where E[ ] is the expected value. The parameter k has been observed for 
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various locations (Eagleson and Wang, 1985) but its choice will be discussed in 
more detail in Chapter III. 

The distribution of surface soil saturation is based on observations by 
Bell et al. (1980) and Owe et al. (1982). As an approximation, the Gamma 
distribution is used: 


f o(s) = 


a 


iR 


o-l -As 
s e 


A,a,s > 0 


where the two parameters are given by 


( 1 . 10 ) 



( 1 . 11 ) 


a — 



( 1 . 12 ) 


where ru g is the coefficient of variation (ratio of standard deviation to the 
mean) of surface soil saturation, and T( ) is the gamma function. Fig. 1.3 
gives field data of cv as a function of area. Note that as area approaches 

o 

c n 

typical GCM grid area (10 km ), cv is approximately 1. 

Surface runoff is produced by two mechanisms: the Dunne mechanism 

and the Horton mechanism. Dunne runoff occurs when the surface soil layer 
is saturated (i.e., s > 1). Horton runoff occurs when the precipitation rate 
exceeds the infiltrability f* of the first soil layer. The infiltrability is given 
by the linearized darcy equation of vertical steady flow in porous media for 
saturated conditions due to ponding: 


f* = K(l) 


dz 


surface saturation 


+ 1 


(1.13) 
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> 


Fijnire 1.3 Dependence of soil moisture variability on element size. 

(Data are from the Hand County, South Dakota site, from Hawley 
et al. (1983), and from the Washita watershed in Oklahoma.) 
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where K(l) is the saturated hydraulic conductivity (i.e., for s = 1), ^ is the 
matric head potential. Application of the chain rule to (1.13) results in 


f* = K(l) 


ds 


s=l 


s -1 
Az 


+ 1 


(1.14) 


Moisture retention curves are used to determine the value of . Here Az 

is the thickness of the surface layer. We are now in a position to write an 
expression for runoff over the wetted fraction k, including both Dunne runoff 
(integrated over all P and s > 1) and Horton runoff (integrated over P > f* 
and 0 < s < 1): 


Q = *[ ff p f p (P) dP f s (s) ds + J J"(p - i*) yp) dP f s ds] (1.15) 
10 Of* 


Substitution of (1.9), (1.10), and (1.14) into (1.15) yields an expression for the 
dimensionless runoff ratio (R = Q/E[P]): 


R = 1 - 



L « l 

7 

“» E 

m 




-«I( l-v) 7 


a, /civ + 


+ 


a 

MiL 


ftlvEfs] 

a 


a 


(1.16) 


+ 1 


T(a) 


where 



1 

Si 


(1.17) 


I = K(1)/E[P] 


(1.18) 


and t(,) is the incomplete gamma function. E[P] is the grid precipitation rate 
(produced by model atmospherics) and K(l) is the saturated soil hydraulic 
conductivity (dependent on soil type). 
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2. EvaTvrtra.nKpiration Efficiency 


As with the GISS parameterization, evaporation is taken as a fraction, /?, 
(called the evaporation efficiency), of potential evaporation (see Eq. (1.5)). 

Spatial variability is introduced into the large grid evaporation estimation 

through soil moisture. 

The derivation of bare soil evaporation efficiency 0 S begins by using the 
desorption equation of Philip (1957), which is an approximation to the basic 
governing equation of vertical moisture diffusion and mass conservation. 


f e “ 2 S e ‘ 


- 1/2 


[K(s ) + K(s,)] 


(1.19) 


where S is a constant based on soil properties and s Q and are initial and 
6 

boundary conditions given by 


s(z,0) = s Q 
s(0,t) = s x 


z,t > 0 


( 1 . 20 ) 


A uniform moisture content s^ holds for the entire profile at relative time 

zero; a constant s ^ content holds at z = 0 at all times. The functional form 

of S is derived by Entekhabi and Eagleson (1989a) based on the work of 
e 

Parlange et al. (1985): 


S 


e 


8nmK( 1)^1) 

3(1 + 3m)(l + 4m) 


1/2 


( 1 . 21 ) 


where n is porosity and m is soil pore disconnectedness index. 
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Unsaturated hydraulic conductivity is parameterized similarly to Brooks and 
Corey (1966): 


K(s) = K(l) s m 

For s 1 < < s Q , (1.19) may be approximated by 

f e * 5 s e ^ + I R ( s o) 


( 1 . 22 ) 


(1.23) 


Under the circumstances of desorption this is a reasonable assumption. 

To accommodate the time step of GCMs (generally one hour), the 
desorption must be integrated to obtain a mean time-step value: 


f 


e 



s e T -1/2 - i K (s„) 


(1.24) 


This equation may be rewritten by substitution of Eqs. (121) and (1.22) to 
yield 


x + 2 . - + 3 

f e = K(1) n s 0 2m - |k( 1) s 0 m 

(1.25) 

o _ \ 8nil>(l) l 1 / 2 

11 ~ [3K(1)T(1 + 3m)(l + 4m) 

(1.26) 


With the mean desorptivity defined, the bare soil evaporation efficiency 
may now be obtained following a derived distribution. Whenever f is less 
than the potential evaporation rate e p , the value of f g is the evaporation rate. 
This is the "soil-controlled" evaporation regime. For f greater than e p , the 
value of e p determines the evaporation rate; this is the "climate-controlled" 
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evaporation regime. A transitional value of relative surface soil saturation s*, 
is defined for which the value of f and ep are equal (i.e., f^ — e p)> 

e s=s* 

From (1.25) we have 

T ~+ 2 ~ +3 

e p = K(l)fi s* 2m - |K(1) s* m . (1.27) 

Using these definitions, then, the expected value of area evaporation or 
the grid evaporation rate is 

s* 

E[e] = e p + f g (s) ds + | f 0 f g (s) ds (1.28) 

s 0 

Finally, integration of (1.28) in combination with (1.25), (1.28), and (1.10) 
results in 



We now move to the derivation of the vegetation transpiration efficiency 

I The derivation is similar to that of (3 except that a soil moisture 
v‘ h 
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extraction function e v (s) by roots takes the place of f . As a simplification of 
a complex system, the following definition of e y (s) is used 


0 


0 < s < s 

W - w 



0 


S < S < S* 

W 

s* < s < 1 
saturated 


(1.32) 


where s is the plant wilting point soil moisture value and s* is that above 
which the transpiration by plants is limited by the climate-controlled potential 
evaporation. The determination of s w is based on the Brooks & Corey (1966) 
equation for partially saturated matric suction: 


H s) = V<1) s 


m 


(1.33) 


This may be rearranged and adapted to the wilting point by 


s = 
w 


t 


wilt 


rm 


-m 


(1.34) 


The derived distribution of evaporation efficiency for vegetated regions 
now follows from (1.32): 


E[e v ] 


e v ( s ) f s (s) ds + E p | f s (s) ds (1.35) 


W 
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Substitution of (1.32) and (1.10) into (1.35) then yields upon integration: 


A, = 


E t e v l 


v e T 


7(a+l,af" 1 )-^" 1 7(a,a£~ i )-7( a + 1 . L )+°M L ^{a,o^ *) 

T(a) (a£ _1 - af -1 ) 


- 1\ . _1U 1—/— «.ll) 


1 + 


(1.36) 


where 


r= 


_ Mi] 


(1.37) 


w 


3 Summary 

Thus the runoff coefficient R and the evapotranspiration efficiency $ 
(with differentiation now between bare-soil and vegetated conditions) have been 
derived from physically based equations of soil water movement in conjunction 
with spatial variability in precipitation and soil moisture. Empirical relations 
have only been invoked in situations such as root extraction of soil moisture 
for which eloquent governing equations simply do not exist. 

Chapter III contains plots comparing these functions to those of the 
current GISS GCM hydrology (see Figs. 3.1a-c). For a more detailed 
discussion and graphical analysis, the reader is directed to Entekhabi and 
Eagleson (1989a). 

D. Screening Model Description and Enhancements 

L Basic Model Description 

The one-dimensional Screening Model used here for further sensitivity 
testing of landsurface hydrological parameterizations has been described in 
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detail elsewhere (Koster et al., 1988, and Entekhabi and Eagleson, 1989b). As 
mentioned earlier, the prime advantage of low-dimension climate models is 
their ability to account for the major feedbacks and sensitivities in the system 
(in this case those of the landsurface-atmosphere) while requiring only a 
fraction of the computer resources of GCMs. 

The Screening Model has one landsurface "grid-square" and one ocean 
surface "grid square" (quotes are to indicate the model does not necessarily 
represent a particular location on earth). Based on the GISS GCM, there are 
nine atmospheric layers over each surface with pressure levels centered at 984, 
917, 807, 650, 480, 329, 206, 105 and 27 millibars. Due to the inherent lack 
of lateral atmospheric dynamics (because of the one-dimensionality), the 
pressure tendency is zero (i.e., dP s /dt = 0). The atmospheric interactions of 
moisture and heat transfer between the two columns are parameterized. After 
Entekhabi and Eagleson (1989b), these transfers follow a simple linear reservoir 
scheme for each atmospheric level i: 


where, per unit mass, 


d Joi 

dt 


f ci - f oi 


df . 

Cl 

dt 


f oi - f d 


(1.38) 



q. 


c T 
C P A i 


moisture exchange 
heat exchange 


The variables q. and T. are the specific humidity and temperature at level i 
and Cp is the specific heat of air at constant pressure. The parameter x is 
in units of time and its effects on moisture and heat convergence in the model 
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have been analyzed in Entekhabi and Eagleson (1989b) 3.nd comp orisons of 
these effects with observations in the tropics have been studied in Brubaker 
et al. (1991). With a dynamically linked system of ocean— atmosphere-land 
there is thus a complete hydrologic cycle operating in the Screening Model. 

The other major parameterizations in the model include rainfall 
production (in the forms of moist-convective origin and large-scale 
supersaturation), radiation, and landsurface hydrology. A brief description of 
each follows here. 

Large-scale condensation occurs by condensing all vapor necessary to 
bring atmospheric saturation down (if need be) to 100%. Supersaturation may 
result from moisture and heat convergence. Moist— convective events are 

simulated by raising subgrid scale moist parcels pseudo— adiabatically from all 
nine atmospheric levels to their points of neutral buoyancy. Any moisture 
which condenses in this process and is not subsequently evaporated at lower 
levels (as droplets fall) becomes surface rainfall. Latent heat transfers are also 
accounted for in this scheme. The parcel fraction used is 1/2 the grid square 
air mass at each atmospheric level. 

The radiation scheme is based on Hoffman (1981). Solar radiation is 
integrated over the model time step and supplies the forcing for the seasonal 
and diurnal cycles in the model. Clouds are simulated for both reflection and 
transmission of radiation. Also the major radiatively active gases are 
represented with their respective absorption spectra. These gases, in order of 
significance, are H 2 0, C0 2 , 0 3 , and 0 2 - 

Sensible heat and evaporation fluxes are parameterized by "drag laws" 
using a prescribed wind speed. The transfer coefficients are functionally 
related to the vertical gradients of temperature and specific humidity. 
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A surface layer (between the ground and the first atmospheric layer) is used 
to give a finer resolution. At each time-step, iterations are performed until 
the heat and moisture fluxes from the ground to the surface layer are equal to 
those between the surface layer and the first atmospheric layer. 

Another facet of the Screening Model which is of importance to the work 
here is the incorporation of the Abramopoulos et al. (1988) multiple soil layer 
moisture diffusion scheme. We present here a summary of the equations 
including a transpiration sink term from each soil layer. 

Each soil layer has thickness Az^ porosity ni, and relative soil saturation 
Si. The stored water depth w. is, 

w. = n.s.Az. (1.39) 

Refer to Fig. 1.4 for a schematic diagram. The fluxes F. entering and leaving 
each soil layer follow the darcy equation 

H — H 

F i = K -r r- 1 (1-40) 

where H is the total hydraulic head, K is the hydraulic conductivity, and the 
subscript i refers to layer i (in the case of H) and the interface at the top of 
layer i (in the case of F and z). Hydraulic conductivity is assumed constant 

in the vertical. The fluxes of evapotranspiration e are evaluated by 

i 

e v = e A e ? ( L41 ) 

i i 

where £i is the root fraction in layer i (the sum of all ei is then unity), /? is 

i 

the evaporation efficiency (evaluated at soil saturation Si) discussed in previous 
sections. Resistances along the plant root system are considered negligible. 
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Figure 1.4 A three-layer version of the Abramopoulos et al. (1988) soil 
moisture diffusion scheme. 
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* 


The change in stored water w- at each soil layer i during a time step At 
is then 


Aw. 

1 

AT 


= F. - F. 


i+1 


V. 

1 


(1.42) 


An implicit time scheme is used for the fluxes resulting in the following 
matrix of linear algebraic equations 


AAw = b 


where 


dF. 


x,i-l dw . 


At 


i -1 


(1.43) 


(1.44a) 



dF. 


i+i 


7W . 
1 


dF. deT.l 



> *J 


(1.44b) 


dF. 


k,i+l 5w. 


i+1 


At 


i+1 


(1.44c) 


b- = (F. - F + e )At 

1 V 1 1-1-1 V/ 

1 


(1.44d) 


The boundary condition at the surface is given by the net moisture entering 
the soil column: 


Fj = P - Q 
= P(i - R) 


(1.45) 


39 



At the lowermost boundary a zero flux boundary condition is used. Sensitivity 
to this will be examined in Chapter II using conditions of an open column 
draining by gravity. 

2 F.nhancPTnents to the Screening Model 

The main modification to the Screening Model has been the incorporation 
of multiple soil layers to simulate seasonal heat storage in the ground. Heat 
storage in the soil is necessary to allow long-duration simulations with a 
seasonal cycle (heretofore the Screening Model has been used mainly for 
perpetual-day simulation). 

To provide seasonal heat storage, we follow Hansen et al. (1983) using a 
finite difference scheme to solve the heat conduction equation. The governing 
equation is 

5 A H =k £T (1.47) 

^ 5z 2 g & 2 

where T is temperature, A is the soil thermal conductivity, is the heat 
capacity per unit volume, and K is the thermal diffusivity. Fig. 1.5 gives a 

o 

soil profile schematic showing the fluxes of heat. For a 3— layer soil scheme, 
the flux equations are 

dT, 

ir = F o - F i (148a) 

*V. ar - F 1 - f 2 (1 ' 48b > 

dT 

ir = f 2 - f 3 < 148c > 


40 






The boundary conditions are given by the net energy balance at the surface 

for F q at each time-step. The lowermost boundary is considered a no-flux 

boundary (F = 0). Using a centered difference scheme, we have for F 1 and 
3 


3T, - 3T„ 


— - F 
2 0 


Az i 

XT 


F, = 


T 


1 F - i 

2 2 A 2 




A z , 


+ 


(1.49a) 


l Az 2 

3T 2 - 3T, - | F, t - 


F 2 = 




Az. 


+ 


(1.49b) 


We assume constant thermal properties for all layers using A - 0.42 w/m/K 
and c = 1.13 * 10^ J/m^/k. 

s 

We have now given the background material necessary for the climate 
simulations done with the Screening Model (Chapter II) and the full 
three-dimensional GISS GCM (Chapter III). Summaries of the results along 
with suggestions and recommendations for the GISS model and for further 
research will be given in Chapter IV, the concluding chapter. 
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Chapter II 


Numerical Experiments with the One-Dimensional 

Screening Model 


A. Introduction 

This chapter is devoted to sensitivity experiments performed with the 
One-dimensional Screening Model described in Chapter I. The effects of 
changes in some key boundary conditions have not been analyzed in the past 
(Entekhabi and Eagleson, 1989b). Among these are the sensitivity to soil 
storage capacity and heat capacity (i.e., depth of the soil column), sensitivity 
to percolation from the lowest soil layer, and sensitivity to transpiration 
extraction from lower soil layers. 

The experiments with the Screening Model are intended to provide 
guidelines for the implementation of the Entekhabi and Eagleson (1989a) and 
the Abramopoulos et al. (1988) hydrologic parameterizations. Sensitivities that 
are determined through experimentation with the Screening Model form the 
basis for the choice of parameters to be used in the full GCM. This approach 
allows more efficient use of the GCMs and their associated computational 
facilities for more focused numerical climate experiments. 

B. Screening Model Experiments 

Here we examine hydrologic sensitivity to soil water and heat capacity, 
to groundwater percolation, and to lower layer transpiration by plants. The 
series of numerical experiments defined for the Screening Model are listed in 
Table 2.1. 


43 



Table 2.1 Simulations performed in the Screening Model 


K— 1 Control Run, light Soil 

K-2 Control Run, heavy soil 

K— 3 Proportionally 1/2 * nominal soil layer thicknesses, light soil 

K _4 Proportionally 2 x nominal soil layer thicknesses, light soil 

K _5 Proportionally 4 * nominal soil layer thicknesses, light soil 

K-6 Top layer: Nominal; lower layers: 4 * nominal thicknesses, light soil 

K— 7 Groundwater percolation from lowest soil layer, light soil 

Groundwater percolation from lowest soil layer, heavy soil 

K_9 Transpiration Control Run, shallow soil column, light soil 
K— 10 Root fractions /.75, .15, .10/ , shaUow soil column, light soil 

K— 11 Root fractions /.85, .10, .05/ , shallow soil column, light soil 

K— 12 Transpiration Control Run, deep soil column, light soil 

K— 13 Root fractions /.75, .15, .10/ , deep soil column, light soil 

K— 14 Root fractions /.85, .10, .05/ ., deep soil column, light soil 

K— 15 Transpiration Control Run, shallow soil column, heavy soil 

K— 16 Root fractions /.75, .15, .10/, shallow soil column, heavy soil 

K— 17 Root fractions /.85, .10, .05/, shallow soil column, heavy soil 

K-18 Transpiration Control Run, deep soil column, heavy soil 
K— 19 Root fractions /.75, .15, .10/, deep soil column, heavy soil 

K— 20 Root fractions /.85, .10, .05/, deep soil column, heavy soil 
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In order to use the hydrological parameterization of Entekhabi and 
Eagleson (1989a), relevant parameters must be assigned representative values of 
their mean over the landsurface region. Table 2.2 contains values for soil 
hydraulic properties as given in Entekhabi and Eagleson (1989a); this 

information is also used in these simulations. In addition, the model requires 
a number of other boundary conditions and specifications. These are given in 
Table 2.3, and are also based on Entekhabi and Eagleson (1989b). 

1. Sensitivity to Soil Water and Heat Capacity 

Table 2.3 designates a "nominal" distribution of soil layer thicknesses for 
a three-layer version of the Abramopoulos et al. (1988) soil moisture diffusion 
scheme (Control Run K-l with "light" soil). These thicknesses follow a 
roughly geometric progression and combine to form a total soil column depth 
of 75 cm. Obviously there are many regions on earth having a mean depth to 

bedrock greater than (or possibly less than) 75 cm. The reason a particular 

model depth is chosen is usually based first of all on the thickness of the 
upper layer since it has a strong impact on fluxes in and out of the soil 

system. Its chosen thickness is also tied to the model time step because often 
an assumption is made that moisture and heat become uniformly distributed in 
the upper soil layer within the period of a time step. Once upper layer 
thickness is determined, a progression is chosen for subsequent underlying soil 
layer thicknesses with geometric progression; this also gives a numerically 
stable behavior. 

To test the effect of soil layer thicknesses, we examine cases wherein the 
proportions of the three soil layers to each other remain the same but the 
total depth is changed to l/2x, 2*, and 4* the nominal total depth. Also, a 
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Table 2.2 


m 

K(l) 

m 

n 

\Ir 

wilt 


The definition of soil hydraulic 
Eagleson, 1989a) 

Pore-size distribution index 

Saturated hydraulic conductivity 
[10 ' 3 m hr" 1 ] 

Saturated matric potential [m] 

effective porosity 

wilting matric potential[10 2 kP a ] 


properties 

(after 

Entekhabi 

SAND 

SILT 

CLAY 

3.3 

1.2 

0.4 

7.5 

2.2 

0.3 

0.23 

0.46 

0.93 

0.25 

0.35 

0.45 

-15 




Table 2.3 Representative specifications for the Screening Model 


PARAMETER 

Simulation Period 
Latitude 

Fixed Ocean Surface temperature 
Cloud temperature standard deviation 

Surface Wind Speed 
Albedos 

Soil texture fractions 
"light" 

"heavy" 

Nominal 3-layer Soil Column 
Thickness 

Rainfall fractional wetting, K 
Soil moisture coefficient 

of spatial variability, CV ' * 
Land-Ocean atmospheric exchange 
parameter, X 


VALUE 


One Year 

(with one 

year spin— up) 

15° North 



25° C 



3° C 



LAND 


OCEAN 

2 ms' 1 


5 ms' 1 

.25 


.35 

SAND 

SILT 

CLAY 

.75 

.20 

.05 

.15 

.15 

.70 

Laver 1 

Laver 2 

Laver 3 

10 cm 

15 cm 

50 cm 

0.60 



1.0 



1.67 days 
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final run using the nominal top layer thickness (10 cm) and 4* nominal lower 
layer thicknesses will be used to test the influence of the top layer thickness 
alone. Light soil textures (see Table 2.3) are used for these simulations. 
Similar results are expected for heavy soil texture since the main sensitivity is 
to heat and moisture storage capacity. 

The heat and moisture diffusion in the soil column is estimated using a 
finite difference scheme with three nodes only. The reasoning behind such a 
coarse discretization is to gain computational efficiency. Deep soil storage of 
heat and moisture is believed to dampen the amplitude of the seasonal cycle 
of modeled hydrologic fluxes as well as the temperature at the surface. The 
total depth of the soil will affect the amplitude of the annual cycle while the 
thickness of the topmost layer should more directly affect the amplitude of the 
diurnal cycles. With the Screening Model, we intend to verify and quantify 
these sensitivities. 

Each of the simulations K-l through K-20 (i.e. , inclusive of the soil 
layer thickness sensitivity simulations considered now) was given a one-year 
spin-up period and diagnostics were collected during the following one-year 
period. Hydrologic equilibrium is reached by the end of each spin-up period 
since, unlike GCMs that have long memory in the ocean component and low 
frequency atmospheric regimes, the Screening Model reaches a stationary state 
once the soil heat and moisture storages are compatible with the net fluxes 
across the land surface. 

We point out that Run K-5 (4* nominal soil layer thicknesses) is 
actually inconsistent with the runoff derivation of Chapter I. This experiment 
is considered to be instructive as an extreme condition. The runoff 

parameterization infiltration equation (Eq. 1.13) requires the top layer 
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thickness to be comparable to the penetration depth of a wetting front over 
one time step (one hour). Run K-5 uses a 40 cm top layer which is too 
thick to meet this condition. 

Finally, before presenting results, we point out that in regard to diffusion 
in the soil, the direct effect of increased soil layer thicknesses will be to 
decrease the fluxes of moisture and heat transfer within a given layer. This 
becomes evident when considering that the moisture flux given by Darcy’s Law 
(see Eq. 1.40) is directly proportional to the total hydraulic head, which is the 
sum of the gravitational and pressure heads: 

H = z + V(s) (2-1) 


Equation (1.40) may then be rewritten 


Fj = K 


- (*i + *,) 


z. , - z. 

i-l i 


= K 


A ^ + Az^ 

Az < 


(2.2) 


where subscript l refers to a layer whose top and bottom elevations are z._ 1 
and z., respectively. Thus since A will be non-zero for virtually all cases 
(except perfectly uniform vertical moisture profiles), any increase in Az^ will 
result in a reduced moisture flux. That is, for a constant nonzero A^, the 
ratio of (A^ + Az^) to (Az^) decreases with increased Az^. By the same 
argument using Eqs. (1.48) and (1.49) the heat flux is also reduced, all other 
variables being constant. 
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The results of the sensitivity simulations for various soil layer thicknesses 
are given in Table 2.4. The diagnostics of the water and heat budgets are 
annual averages and ranges (in parentheses) of: 


8 ( 1 ) 

P 


e 

Q 

-V-q 

T 

g l 

LE 

SH 

F 

g 

-V-F 


Relative soil saturation of top layer [ ] 

Precipitation [mm/day] 

Mean precipitation intensity (i.e., intensity when 
precipitating) [mm/day] 

Evaporation [mm/day] 

Runoff [mm/day] 

Atmospheric convergence of moisture 
over landsurface [mm/day] 

Temperature of top soil layer [C] 

Latent heat flux over landsurface [W/m 2 ] 

Sensible heat flux over landsurface [W/m 2 ] 

Heat flux into soil [W/m 2 ] 

Atmospheric convergence of heat 
over land surface [W/m 2 ] 


Two salient features of Table 2.4 are that 1) as the total soil column 
depth is increased and the layer thickness proportionality is maintained the 
climate becomes steadily and considerably drier (compare s(l), P, E and R as 
well as temperatures in Runs K— 1 to K— 5), and 2) When the thickness of the 
topmost layer is held constant, even large changes in lower layer thicknesses 
produce little or no change in the mean annual climate (Run K-6 and 
Run K-l are nearly identical); the annual range of climatic variables, however, 
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K-4 (2* Nominal) 25.3 (6.8) 89.3 (42.5) 29.1 (42.4) -0.4 (3.9) -12.0 (58.8) 

K-5 (4x Nominal) 27.2 (10.3) 78.7 (36.7) 29.9 (45.1) 0.4 (3.8) -12.6 (72.8) 

K_6 ( lowei- N Layers: 4x> 23 3 (*•*) 105 -° ( 32 ‘ 6 ) 23 8 ( 291 ) 0 0 ( 3 3 ) *«- 4 (58-7) 



is generally decreased as the lower layer thicknesses are increased. Supporting 
the first point we see that as the total soil depth increases from 1/2* nominal 
to 4* nominal depth (Runs K-3, K-l, K—4, K— 5), all components of the 
annual mean water balance are decreased. In the heat balance we see a steady 
increase in mean surface ground temperature and a steady decrease in latent 
heat flux. As drier conditions ensue the sensible heat flux takes on a greater 
role in the surface cooling since there is less water available for cooling by 
latent heat flux. As a measure of aridity, the Bowen ratio of sensible to 

latent heat flux generally follows a rise with greater total soil depth. 

There is a striking similarity of Run K-6 to Run K-l in all of the 
annual mean diagnostics (Run K-6 has a nominal top layer depth of 10 cm, 
with lower layers 4* the nominal depth). The annual range , however, 
particularly evident in s(l) and T ., is reduced by the larger storage capacity 
in Run K-6. One other obvious feature of Table 2.4, to be discussed shortly, 
is the dramatic increase in precipitation intensity P- when the climate becomes 
warmer and drier. 

It was mentioned earlier that the total depth of the soil ought to affect 
the annual cycle of heat and moisture while the thickness of the top layer 
ought to have a greater impact on the diurnal cycle. Thus Run K-6 ought to 
follow the trend of the other simulations in annual range according to its total 
soil column depth while following the diurnal range according to its top soil 
layer depth. This is borne out in Fig. 2.1, which demonstrates this hypothesis 
in graphical form with respect to temperature by plotting the mean annual 
range of surface ground temperature against total soil column thickness (top) 
and the mean diurnal range against the top soil layer thickness (bottom). As 
shown, Run K-6 follows the trend of the others in both plots. 
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Figure 2.1 
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There is an apparent discrepancy in the upper plot of Fig. 2.1 in the 
case of Run K-5. In spite of having the highest heat and moisture capacity, 
it has a very large mean annual surface ground temperature range. The 
reason for this lies in the secondary effects (feedbacks) which are producing 
excessively dry and hot summer conditions. In the winter, the soil moisture 
increases again to give lower surface temperatures roughly comparable to those 
of Run K-4. Thus the annual range of Run K-5 is quite large. 

Annual ranges of relative surface soil moisture are given in Fig. 2.2. 
Again, for the annual range of s(l) we see Run K-6 behaving according to its 
total soil column depth. Recall from Table 2.4 that the annual mean of s(l) 
in Run K-6 is nearly identical to that of Run K-l (control). This is also 


true of annual mean surface ground temperature T 


gl‘ 


Thus we may conclude that the climate mean values of the water and 


heat balance diagnostics are influenced mainly by the top soil layer thickness. 


The diurnal ranges also depend most strongly on the top layer thickness. The 


annual ranges, however, depend on the total soil column depth. 


2. Moist Convection Intensity and Land Surface Temperature 

As evident in Table 2.4, as precipitation in the mean is reduced and the 
climate becomes warmer and drier, precipitation intensity is nonetheless greatly 
increased. By precipitation intensity we mean the average rainfall rate when 
it is raining. Rainfall in the Screening Model is produced mainly by the 
moist-convective mechanism, which is dependent on vertical potential 
temperature and specific humidity gradients in the atmosphere. Fig. 2.3 shows 
that the increase in mean precipitation intensity P. is strongly dependent on 
mean surface ground temperature. Moist static stability, the main criterion of 
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Precipitation Intensity (mm/day) 


Precipitation Intensity vs. 
Surface Ground Temperature 
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Surface Ground Temperature (C) 


Figure 2.3 Mean annual precipitation intensity versus surface ground 
temperature 
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the moist-convective mechanism, is apparently in direct correlation with 
surface ground temperature for these ranges since the ground temperature is an 
indicator of energy influx at the bottom boundary of the air column. This 
energy drives the atmospheric overturning and convection. 

3. Sensitivity to Groundwater Percolation 

This section examines sensitivity to a simple formulation of groundwater 
percolation from the soil column as given by Runs K-7 and K-8; Runs K-l 
and K-2 serve as control cases. 

In most GCM landsurface hydrological parameterizations a no-flux 
boundary condition is used at the bottom of the soil column for the solution 
of soil moisture diffusion equations (see Fig. 1.3). It may be advantageous 
however, particularly in humid regions or places with deep soils, to modify the 
standard soil column to account for groundwater runoff. As a test of climate 
sensitivity to groundwater runoff we have simply changed the lower boundary 
condition of Fig. 1.3 (using Eqs. 1.42 through 1.44d) from a no-flux condition 
to one of a constant flux through the lower layer by setting F 4 equal to F 3 - 
Percolation through the lower boundary is then considered groundwater runoff. 

Two simulations were run--one using light soil (Run K— 7) and one using 
heavy soil (Run K-8). All other boundary conditions and specifications are 
given by Tables 2.2 and 2.3. A 1-year spin-up period was used followed by a 

1-year period of diagnostics collection. 

The results are presented in Table 2.5, using the same notation as that 
for Table 2.4 except that Q T0T represents the total runoff and Q GW 
represents the groundwater runoff. These runoff diagnostics show that the 
strongest impact of groundwater percolation is in the light soil case. Thi? is 
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Table 2.5 Annual mean water and heat budgets for simulations 
showing sensitivity to groundwater percolation 
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Heavy Soil 

K-2 (control) 23.1 (5.8) 104.3 (35.5) 22.7 (26.6) -0.3 (2.9) -7.7 (44.2) 

K-8 (percolation) 23.3 (5.8) 109.7 (39.3) 20.9 (15.9) -0.3 (2.7) -10.3 (56.2) 



reasonable since light soil has a much higher saturated hydraulic conductivity 
than heavy soil (see Table 2.2). Although Q qw in the light soil case 
increased from zero (control) to 0.38 mm/day in Run K— 7, nevertheless the 
surface runoff was correspondingly decreased such that total runoff Q'pQ'p 
increased by only 10%. The surface relative soil saturation s(l) was 
nonetheless decreased substantially by (roughly 30%) since a greater portion of 
soil moisture was partitioned to the lower layers. Because of the drier top 
layer, the surface fluxes of precipitation, evaporation and runoff all decreased 
from those of the Control Run. The reason the total runoff Q'j'Q'p increased is 
that the groundwater runoff more than makes up the difference for the reduced 

surface runoff. 

In the heavy soil case (Run K-8), groundwater runoff is practically 
negligible. It does, however, have a slight effect on the modeled climate. 
Precipitation and evaporation are both increased over the Control Run 
although the surface runoff remains unchanged and the mean soil moisture 
actually decreases. The annual ranges of the water balance variables, 
particularly for precipitation, are all increased. Apparently the mean values of 
P and E are increased due to their high values during humid time periods. 

It is not clear what mechanism is responsible for the marked increase in 
the annual range of P for both the light soil and heavy soil cases. A possible 
explanation is that it is due to the interaction of precipitation and evaporation 
with the soil. Since the top soil layer has less stability with regard to 
moisture retention due to groundwater percolation, it may be that the feedback 
of the atmosphere involves a high variability in precipitation. However, the 
ranges given are of monthly averages so that increases in the range must 
involve long time constants in the system. Also, atmosphere— soil feedbacks 
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ought to affect evaporation in the same way as precipitation, yet we see the 
increase in the range of evaporation is much smaller than the increase in the 
precipitation range. 

We conclude that groundwater runoff in this simple formulation is most 
significant in light textured soils having the effect of increasing total runoff 
while creating an overall drier climate. 

It should be noted that for actual implementation in a GCM, 
topographic information ought to be used to obtain an average terrain slope to 
be used in the groundwater runoff formulation. Groundwater percolation is the 
flow rate under gravity head. For a bedrock with a certain slope, the 

groundwater percolation equals the unsaturated hydraulic conductivity 

multiplied by this slope. Since slopes are generally small over large regions, 

this contribution to large-scale or continental runoff may be small although it 
may play an important role in redistributing moisture in localized parts of 

river basins. 

The experiments performed through Runs K-7 and K-8 represent open 
drainage at the bottom of the soil column (i.e., gravity hydraulic head is 

unity). They thus represent the extreme case and the upper bound on the 
influence of groundwater percolation in GCM landsurface hydrology. 

4. Sensitivity to Lower Laver Transpiration 

This section examines sensitivity to the root distribution in the soil 
column. The experiments are defined as Runs K-9 through K-20 in Table 2.1. 

In the simulations discussed thus far, transpiration has been allowed from 
the top soil layer only. This corresponds to a root distribution of e l = 1.0, 
c 2 = 0.0, e 3 = 0.0 (see Eq. 1.41). As a compact way of representing this we 
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use the notation /1.0, 0.0, 0.0/ for the three soil layers. Because of the wide 
diversity of vegetation over areas as large as GCM grids, exact representation 
of plant root distribution is impossible. Nevertheless, moisture extraction by 
roots can be a significant component of evapotranspiration. To identify model 
sensitivity to these effects a set of simulations is performed as a guide for 
implementing the Abramopoulos et al. (1988) parameterization into the 

GISS GCM (to follow in Chapter III). 

Based on the studies of plant root geometry by Epstein (1973), the root 
distributions used in the simulations are roughly exponential, with the greatest 
concentration of roots in the upper soil layer. Because preliminary simulations 
showed a strong tendency toward drying in the lower soil layers, the fractions 
of roots in the lower layers were kept relatively low. In the simulations 
presented here, we use root distributions of /0.85, 0.10, 0.05/ and /0.75, 0.15, 
0.10/ for cases of heavy soil and light soil. Furthermore, based on the results 
of II.B.l and anticipating the implementation of the soil diffusivity scheme 
into the GISS GCM, two cases of soil layer thicknesses are used— a shallow 
case and a deep case. Both of these cases use a top layer thickness of 10 cm 
and a middle layer thickness of 15 cm. The lowest layer thickness is set to 
25 cm for the shallow case and 200 cm for the deep case. All other boundary 
conditions and specifications are as given in Tables 2.2 and 2.3. 

The results of the simulations are presented in Table 2.6 following the 
same format as Tables 2.4 and 2.5. In general, the effect of increasing the 
transpiration flux from the lower soil layers is a slightly more humid climate. 
The top layer is slightly more moist since the lower layers are now taking 
some of the burden for seasonally high evapotranspiration. While the 
sensitivity of the top soil layer moisture ranges from being essentially 
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Table 2.6 Annual mean water and heat budgets for simulations testing sensitivity 

to lower layer transpiration 

(Numbers in parentheses show annual ranges of monthly averages; 

Numbers following simulation names give 3-layer root distributions) 


Light soil, shallow column 
K-9 (/1,0.,0./) 

K-10 (/.75,.15,.10/) 

K-ll (/.85,.10,.05/) 

Light Soil, Deep Column 
K-12 (/l.,0.,l./) 

K-13 (/.75, .15..10/) 

K-14 (/.85,.10,.05/) 

Heavy Soil. Shallow Column 
K-15 (/l.,0.,0-/) 

K-16 (/.75,.15,.10/) 

K-17 (/.85,.10,.05/) 

Heavy Soil. Deep Column 
K-18 (/l.,0.,0./) 

K-19 (/.75,.15,.10/) 

K-20 (/.85,.10,.05/) 


sGQIl 

Pfmm/dav] 

.17 (.15) 

3.82 (1.85) 

.20 (.15) 

4.09 (1.73) 

.20 (.09) 

4.15 (1.51) 


.15 (.12) 

3.98 (1.67) 

.22 (.08) 

4.47 (1.81) 

.19 (.12) 

4.11 (1.97) 


.20 (.19) 

4.41 (1.86) 

.21 (.18) 

4.45 (2.20) 

.20 (.13) 

4.19 (0.99) 


.19 (.17) 

4.18 (1.49) 

.27 (.11) 

4.69 (2.19) 

.24 (.17) 

4.38 (1.84) 


Elmm/davl 

Qfmm/dayl 

3.28 (1.21) 

0.57 (0.66) 

3.46 (1.28) 

0.64 (0.82) 

3.59 (1.18) 

0.63 (0.38) 

3.14 (0.94) 

0.49 (0.62) 

3.84 (1.79) 

0.89 (0.80) 

3.56 (1.30) 

0.62 (0.55) 

3.67 (1.57) 

0.74 (0.74) 

3.64 (1.42) 

0.80 (0.88) 

3.60 (1.31) 

0.62 (0.55) 

3.58 (1.21) 

0.61 (0.74) 


3.97 (1.90) .08 (0.84) 

3.79 (1.41) 0.85 (0.70) 
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unaffected (Runs K-15 through K-17) to being rather strongly affected (Runs 
K-18 through K-20), the changes in the water balance variables are generally 
within 10% of the control. 

The stronger effect is occurring in the lower layers themselves. 
Figs. 2.4a and b show daily values of the relative soil saturations of all three 
layers for the cases of light and heavy soil in a shallow soil column. As can 
be seen, particularly in the light soil case (Fig. 2.4a), there is a substantial 
drying of the middle and lowest soil layers (i.e. , layers 2 and 3) when 
transpiration from those layers is allowed; It is interesting to note that the 
relative saturation of layer 1, however, is relatively unaffected by the strong 
changes in layers 2 and 3. The main reason for this is that while there is 
less of a water supply from the lower layers, there is also a reduced 
evaporative draw from layer 1 since it has a reduced root fraction from the 
Control Run. 

Based on these simulations, and in the absence of root distribution data 
for the GISS GCM, the choice was made to use a root distribution of 
/.85, .10, .05/ for implementation into the GISS GCM (described in 

Chapter III). The severe drying of lower soil layers found with larger 
fractions of roots in lower soil layers is regarded as unrealistic for most 

settings. 
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Light Soil; Shallow Column; 

Root Fractions /1.0, 0.0, 0.0/ (Control Run K-9) 



Uyer 1 

— Layer 2 
— — Layer 3 


. ?« 

• , v v M- u . ,v . 


Year 1 


Year 2 


Light Soil; Shallow Column; 



Light Soil; Shallow Column; 

Root Fractions /0.75, 0.15, 0.10/ (Run K-10) 
1.0 1 _ 



Year 1 


Year 2 


Figure 2.4a Daily 3-layer relative soil saturations for the case of light soil 
and shallow soil column for root distributions of /1.0, 0.0, 0.0/ 
(top), /0.85, 0.10, 0.05/ (middle), and /0.75, 0.15, 0.10/ (bottom) 
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Heavy Soil; Shallow Column; 



Year 1 Year 2 


Heavy Soil; Shallow Column; 



Year 1 Year 2 


Figure 2.4b Daily 3-layer relative soil saturations for the case of heavy soil 
and shallow soil column for root distributions of /1.0, 0.0, 0.0/ 
(top), /0. 85, 0.10, 0.05/ (middle), and /0.75, 0.15, 0.10/ (bottom) 
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C. Discussion 

This chapter has essentially been a set of screening simulations in 
preparation for the implementation of the spatial variability hydrological 
parameterization of Entekhabi and Eagleson (1989a) and the soil moisture 
diffusion scheme of Abramopoulos et al. (1988) into the GISS GCM. Based 
on the results of Part II.B.l the upper two soil layers will be kept constant 
for all grids with layer 1 thickness 10 cm and layer 2 thickness 15 cm. The 
lowest layer thickness will then be adjusted to give a consistent total soil 
storage capacity with that used in the current GISS GCM. Layer 3 thickness 
will be set such that the total field capacity of the new soil column, based on 
its porosity, will be equal to the field capacity in GISS Model II. 

The results of the groundwater percolation simulations in part II. B. 2 
showed significant runoff sensitivity only for cases of light soil where the 
hydraulic conductivities are generally high. In the absence of a means for 
obtaining representative terrain slope information, the simulations in the 
GISS GCM will use the no— flux lower boundary condition. 

Finally, due to the lack of reliable plant root distribution observations, it 
has been determined that a distribution of /.85, .10, .05/ yields realistic results 
and it will be used in the GISS GCM experiments of the following Chapter. 
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Chapter III. 

Numerical Experiments in the GISS GCM 
A. Introduction 

1. Overview of th* Numeric Experiments 

The goal of the effort described thus far is an improved hydrologic 
scheme whose utility may be demonstrated in a three-dimensional GCM. As 
explained in Chapter I, the current empirical relations in the GISS GCM 
hydrologic parameterization are not physically-based. An area of improvement 
for this GCM would be to introduce more realistic expressions for land surface 
hydrologic fluxes such as those proposed by Entekhabi and Eagleson (1989a). 
Their scheme introduces subgrid-scale spatial variability for some key 
parameters and derives physically-based expressions for the grid-average rates 
of r un off and evapotranspiration. Furthermore, the improved GCM land 
surface hydrology scheme would also require a more realistic soil moisture 
storage and diffusion scheme in order to realistically account for seasonal cycles 
in surface water and heat balance. For this purpose we include the soil 
moisture diffusion scheme of Abramopoulos et al. (1988) alongside the 
landsurface parameterization of Entekhabi and Eagleson (1989a). 

The implementation of the new schemes into the GISS GCM is 
performed in a step-wise fashion in order to isolate the marginal changes due 
to each scheme. Four simulations with the three-dimensional GISS GCM will 
be compared as listed in Table 3.1. The numerical experiments will use the 
8 * io degree GISS GCM with fixed sea-surface temperatures based on 
climate values. Run G-0 is the control case and its diagnostics are collected 
for a three-year period. Run G-l differs from G-0 only due to the 
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Table 3.1 Simulations performed in the GISS GCM. 


Simulation Name Description 

G— 0 GISS Model II as described in Hansen et al. 

(1983) using 8 * 10 degree grid resolution and 
fixed ocean temperatures (Control Run). 

G— 1 Same as G-0, except for new formulation of 

runoff coefficient R, bare soil, evaporation 
efficiencies, /? s , and transpiration efficiency, (3 V 
(after Entekhabi and Eagleson, 1989a). 
Fractional wetting parameter k set equal to 1.0 
for large scale supers at uration rainfall and 0.6 
for moist-convective rainfall. 

(Abbreviated M space M ) 

G— 2 Same as G— 1, except new soil moisture 

diffusion scheme of Abramopoulos et al. (1988) 
is used with transpiration from lower layers. 

No instantaneous upward diffusion or pre- 
scribed growing season as in GISS II. 

(Abbreviated "space/soil n ) 

G-3 Same as G-2, except fractional wetting 

parameter k for moist-convective rainfall set to 
0.15 

(Abbreviated "space / soil/storm" ) 
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replacement of the Hansen et al. (1983) empirical runoff and evapotranspiration 
efficiency functions with the Entekhabi and Eagleson (1989a) expressions. Run 
G-2 is further modified to include the Abramopoulos et al. (1988) soil 
moisture diffusion scheme. The final simulation experiment G-3 differs from 
G-2 only by changes in one parameter that defines the spatial fractional 
wetting in the subgrid parameterization for runoff generation. 

The diagnostics collected for each run in addition to those standard for 

the GISS-II model include: 

a) For all grids; monthly values of potential evaporation, evaporation 
over earth, runoff produced by snowmelt, no-rain probability, 
percentage of rainfall which is of moist-convective origin, number of 
independent storm events, mean storm duration, mean interstorm 
duration, and mean storm depth, 

b) For selected regions; daily values of precipitation, evaporation, 
surface temperature, and relative surface soil saturation, 

and c) For three selected United States grids; hourly precipitation records. 

Initial conditions for each of the simulations are taken from the standard 

GISS Model II 8" * 10" Run as described in Hansen et al. (1983). 
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2. Model Boundary Conditions of the Landsurface 

The GISS-II GCM landsurface boundary conditions are based on the 
archived 1° * 1" vegetation data set of Matthews (1983) and the archived 
1° * 1° soil texture data set of Zobler and Cary (1984). Eight vegetation 
types were defined by Matthews, and each vegetation type was associated with 
a value for albedo, roughness length, "masking depth" (a variable describing 
the effect of vegetation on snow albedo), and field capacity. The values of 
the parameters used over a particular grid in the model are weighted averages 
based on the percentage of each vegetation type over that grid. Table 3.2 
gives the values of parameters for each vegetation type. 

Soil textures (given by percentage of sand, silt, and clay) are assumed to 
have no variation in the vertical dimension. From the soil textures, values of 
the soil hydraulic properties are taken as weighted averages over the three soil 
types. We follow Entekhabi and Eagleson (1989b) in assigning values as 
presented in Table 2.1. 

3. Data Sets Used to Evaluate Model Performance 

Simulated model climates for each of the cases (G-0 through G-3) are 
compared with each other and with some observed hydroclimatic variables. 
The sources and data quality for each of the observed data sets used in the 
comparisons are described below. 

Comparisons will first be made of the water balance over the entire 
globe (divided into lumped values over the ocean and landsurface). From 
there, we move progressively to finer scales of resolution including zonally 
averaged (i.e., latitude-belt averages) values, continentally averaged values, and 
finally balances over major hydrologic basins. Data sets used are those of 
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► 




Desert* 

Tundra 

Grass 

Visual albedo 

Winter 

0.35 

0.07 

0.09 


Spring 

035 

0.06 

0.10 


Summer 

035 

0.08 

0.09 


Autumn 

035 

0.08 

0.09 

Near-IR albedo 

Winter 

035 

0.20 

0.27 


Spring 

035 

0.21 

035 


Summer 

035 

030 

036 


Autumn 

035 

0.25 

031 

Field Capacity 

Layer 1 

10 

30 

30 

(g/tn) 

Layer 2 

10 

200 

200 

Masking Depth 

(m) 


0.1 

03 

0.2 

Roughness Length 

(m) 

0.005 

0.01 

0.01 

Wilting Point 
(10 2 kPa) 


— 

-30 

-22 


Shrub 

Woodland 

Dediduous 

Evergreen 

Rainforest 

0.09 

0.08 

0.10 

0.07 

0.06 

0.10 

0.07 

0.05 

0.07 

0.06 

0.14 

0.08 

0.06 

0.08 

0.06 

0.11 

0.06 

0.05 

0.06 

0.06 

0.27 

0.23 

030 

0.20 

0.18 

030 

0.24 

0.22 

0.20 

0.18 

0.42 

030 

0.29 

0.25 

0.18 

033 

0.20 

0.22 

0.18 

0.18 

30 

30 

30 

30 

200 

300 

300 

450 

450 

450 

03 

2 

5 

10 

25 

0.018 

032 

1 

1 

2 

-20 

-20 

-20 

-20 

-15 


* Desert albedo is reduced by a factor of (1 - 0.5 Si) 


Table 3.2 


Vegetation boundary conditions of GISS Model II 

All values derived from Hansen et al. (1983) except for the 

wilting level values which are based on Entekhabi and Eagleson 

(1989a) 
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Budyko (1978) for global hydrologic balance, Zubenok (1970) for zonally 
averaged balance over landsurfaces, Henning (1990) for balance over continents, 
and Russell and Miller (1990) for major river basin water balances. 

Further comparisons include the seasonal distribution of surface air 
temperature and precipitation over the earth. January and July mean surface 
air temperature data from Washington and Meehl (1984) are used, as well as 
December-January-February (DJF) and June-July-August (JJA) mean 
precipitation data from Schutz and Gates (1971). 

Finally, we consider precipitation regime comparisons for gridded data 
from the National Center for Atmospheric Research (Sprangler and Jenne, 
1988), and rainfall frequency comparisons based on hourly station data over 
the U.S. from the National Climate Center data set. As a check on the heat 
balance, we also compare zonally averaged temperatures over landsurfaces to 
the data of Legates and Willmott (1990). 

Whereas comparisons of the GISS Model climate against observed data 
have in the past been done over the entire globe, the comparisons performed 
here will concentrate on the landsurface only. Comparisons of this sort are 
vital to the validation of parameterizations such as those generating runoff and 
evaporation over the landsurface since in global comparisons that include land 
and ocean, the behavior of the oceans can obscure the dynamics of the 
landsurface. A case example is evaporation; the zonally-averaged evaporation 
rate is dominated by the ocean evaporation rate at many latitudes. When the 
GCM is characterized by seasonally-fixed ocean temperatures, this diagnostic is 
a poor indicator of model performance, especially over land grids. 

Data sets whose goal is to provide global coverage often have poor 
quality in regions where observation stations are sparse. Also, in the case of 
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zonally averaged data over landsurfaces, interpretations must be made keeping 
in mind the percentage of a given latitude belt which is in fact landsurface 
and not ocean. Where this percentage is small (such as in the 40° to 60° S 
Latitude region in Figs. 3.18a-c) observed and simulated data must be viewed 
with greater caution. Throughout the comparisons we attempt to point out 
data quality issues such as these. 

B. Implementation of the New Hydrological 
Parameterization into the GISS GCM: 

Some Off-Line Results 

1 Description 

Incorporation of the new hydrology schemes into the GISS GCM requires 
adaptation to some new details not present in the screening model used in 
Chapter II. In this section, the manner in which the Entekhabi and Eagleson 
(1989a) and Abramopoulos et al. (1988) schemes were incorporated into the 
GISS GCM are outlined. A few off-line (outside of the GCM) sensitivities are 
also analyzed by assuming typical forcing for each of the schemes. These 
off-line tests are important indicators of expected model sensitivities. 

a. Runoff Ratio and Evapotranspiration Efficiency . Precipitation is 
generated in two forms over the grid square: moist convective rainfall and 

large-scale super saturation rainfall. We have taken their occurrences to be 
spatially independent for any given time step in the model. Given an amount 
of each type of rainfall, we use Equation (1.16) to obtain a runoff coefficient 
for the moist convective event, R mc , an d for the supersaturation event, R gg . 
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With precipitation intensities P and P , we obtain the total runoff Q as 

me So 


Q = MAX [P mc R mc + P SS R SS , 

soil moisture capacity exceedance] (3.1) 

The new parameterization must also accommodate the effects of freezing 
temperatures. Specifically, snowmelt runoff and ice in the soil column need to 
be addressed. The current GISS parameterization treats snowmelt runoff 
identically to rainfall runoff. We choose to deal with snowmelt runoff in two 
ways. In Run G-l we allow no runoff from snowmelt except when the field 
capacity is exceeded. In Runs G-2 and G-3 we treat it in the same way as 
supersaturation rainfall, i.e., snow coverage is over the entire grid. For the 
case of G-l, we see runoff volumes in high latitudes generally peak higher in 
spring thaw than in the control case. This is because any snowmelt during 
the winter remains stored in the ground until the capacity is filled. 

Ice in the soil column is treated as follows: it contributes to the 

relative soil saturation used to formulate both the runoff coefficient and 

evaporation efficiency, but is not allowed to diffuse in the soil. Sublimation 

from the surface layer is allowed in this parameterization. 

The GISS GCM has eight vegetation-type designations, and each grid 
has eight values delineating the fraction of each type within the grid. As we 
have mentioned earlier, Table 3.2 gives values of parameters associated with 
each vegetation type (all values except wilting point values are from Hansen et 
al., 1983). The evaporation efficiency is calculated as a weighted average of 
that obtained over bare soil and vegetated fractions, the vegetated fraction 

evapotranspiration efficiency itself being a weighted average taken using all 
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values of the wilting point. Soil hydraulic properties for the sand, silt and 
day classes correspond to those of Table 2.1. 

The fractional wetting parameter k (representing the fraction of the grid 
square actually wetted by a precipitation event) exerts a strong influence in 
the spatial variability parameterization of runoff as seen in Figs. 3.1a-e. 
Unfortunately there is no simple method for the determination of a proper 
value of k for use in GCMs. In nature, k would have a complicated seasonal 
and geographical dependence on regional climate and may be subject to various 
scales of spatial and temporal variability. To provide reasonable estimates for 
k we reduce the complexity of the conceptual model to one simply involving 
storms over grid squares. 

Within an individual grid square, there are at least three scales of spatial 
variation to consider. First, there is the size of a typical storm area relative 
to GCM grid area. For large scale supersaturation events these two may be 
similar; moist-convective storms, however, would generally be much smaller in 
area than the entire grid square. Second, there is the fraction within the 
storm itself which is actually wetted by rainfall (e.g., there are generally 
rain— free gaps in meso— scale storm bands). Observational analysis of air mass 
thunderstorms over Arizona and the Sudan (Eagleson, 1987; Eagleson, 1984) 
suggests this value may be roughly 0.60 in the mean for moist-convective 
events. Third, there is the relative location of storms within the grid square. 
To reflect the spatial variability of storm location, storms may be modeled as 
only partially covering a grid square (see Eagleson and Wang, 1985). The 
value of «, then, ought to be determined based on factors related to these 
interactions. 
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Figure 3.1c Runoff coefficient R for k = 0.1, cv = 1.0, sand (top) and clay 
(bottom) for typical ranges of soil moisture and precipitation 
intensity 

Solid lines = R (spatial variability parameterization) 

Dotted lines = R (GISS parameterization) 






Potential Evaporation (mm/day) 


Figure 3.1d Bare soil evaporation efficiency /3 S for cv = 1.0, sand (top) and 
clay (bottom) for a typical range of potential evaporation 

Solid lines = /3 s (spatial variability parameterization) 

Dotted lines = 3 S (GISS parameterization) 






Figure 3.1e Vegetated soil transpiration efficiency /3 V for cv = 1.0, sand (top) 
and clay (bottom) for a typical range of potential evaporation. 
The wilting pressure is taken as ip = -15 bars 

Solid lines = /? v (spatial variability parameterization) 

Dotted lines = /3 V (GISS parameterization) 


As a simple approximation, we have prescribed one value of k for 
supersaturation and one value for moist-convective rainfall. Due to the large 
spatial scales associated with supersaturation events, a value of k = 1.0 is 
used for this type of precipitation for all simulations. For moist-convective 
rainfall, Runs G-l and G-2 follow Entekhabi and Eagleson (1989a) using 
k = 0.60. In light of the previous discussion, however, it is possible that k 
ought to be set at a much lower value. To test this sensitivity, « for 
moist— convective rainfall is set to 0.15 in Run G— 3. 

One final point is that for Kuo-type moist-convection parameterizations 
(not used in GISS Model II), the fraction of the grid square experiencing 
convection is a variable that depends on boundary layer properties. The 

landsurface hydrologic parameter « may thus be tied to and coupled with the 
moist-convective parameterization if the Kuo-scheme is used in the GCM. 

The contour plots of runoff ratio and evaporation efficiency reveal major 
differences between the linear relationships in the current GISS-II formulation 
(which has dependence on relative soil moisture alone) and the new spatial 
variability formulation (which has nonlinear dependence on relative soil 
saturation, precipitation intensity [in the case of runoff coefficient] and 
potential evaporation [in the case of evaporation efficiency]). 

The nonlinear dependence of the runoff coefficient on precipitation 
intensity is of importance. In Figs. 3.1a-c for a given soil saturation, the new 
parameterization generally gives a lower value of R for low intensity rainfall 
than the GISS formulation. However, as precipitation intensity increases at 
this soil saturation, the new formulation for R exceeds the GISS value. This 
nonlinear dependence is consistent with equations of infiltrability and therefore 
has greater physical validity than does the GISS parameterization. 
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The rather strong dependence of R on k can be seen by comparing 
Figs 3.1a through 3.1c. The runoff ratio for the new parameterization is 
highly sensitive to the parameter representing fractional wetting over the grid. 
The effects are strongest mainly for higher values of precipitation intensity. 
As we will see, the bulk of the grid squares have low-intensity precipitation 
and the effect of the spatial variability will be small. However, many of the 
grids with high-intensity precipitation are tropical grids for which (the Amazon 
River Basin, for instance) current model runoff values are known to be low. 
Fig. 3.2 shows the percentage of the annual rainfall in the Control Run G-0 
that is due to moist convection. As expected, the tropical regions have the 
highest percentage values (up to 70%) and they are hence characterized by 
greater rainfall intensities. In these regions the changes resulting from the new 
runoff parameterizations will be clearly evident and significant changes in 
model climate are to be expected. In these same regions, the sensitivity to 
the landsurface parameter k is going to be highest as evident in Figs. 3.1a 
through 3.1c. 

The evapotranspiration efficiency functions /? g and 0 are plotted in 
Figs. 3. Id and e respectively as functions of the relative soil saturation and 
the potential evaporation rate. These functions are most sensitive to the 
potential evaporation rate at low values of this variable and have reduced 
sensitivity at high values. 

b. New Three-Laver Soil Column (Runs G-2 and G-3) The two-layer 
soil moisture diffusion scheme of the GISS GCM has major shortcomings. 
First, upward diffusion from layer 2 to layer 1 is restricted to a pre-defined 
growing season in mid and upper latitudes (Hansen et al., 1983). The second 
layer is provided to overcome excessive drying of layer 1 in continental regions 
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CONTROL RUN (G-0) % MOIST-CONVECTIVE PPT 



Figure 3.2. Percentage of rainfall which is of moist-convective 
origin in the GISS Model II Control Run 
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(luring summer. The lower layer can recharge outside the growing season and 
then supply its store of water to layer 1 over the summer. The simplified 
differential equations governing diffusion in the current GISS GCM are 
themselves not physically based. For example, there is an instantaneous 
equilibration of relative soil moisture in layers 1 and 2 for vegetated fractions 
of the grid when diffusion is in the upward direction. As an improvement on 
this, we use a three-layer soil column based on the work of Abramopoulos et 
al. (1988). 

This scheme is a finite difference solution of the diffusion equation; we 
use a no-flux boundary condition at the base of the column. A flux boundary 
condition at the top is given by the net inflow into the soil resulting from 
precipitation, evaporation, and runoff. Transpiration from lower layers is 
performed by means of a sink term in each layer proportioned by the assumed 
fraction of roots in each layer. A schematic comparison of the diffusion 
scheme in the current GISS GCM and the one implemented here is given in 
Fig. 3.3. 

The determination of soil layer thicknesses was in part based on the field 
capacities used in the GISS model. Porosities based on the weighted soil 
textures were used to scale the total depth of the column such that the total 
field capacity remained unchanged from the GISS-II GCM. However, due to 
the strong influence of the thicknesses of the upper layers on infiltration and 
desorption rates, the upper two layer thicknesses were set uniformly over the 
entire earth. The top layer is 10 cm and the middle layer is 15 cm. The 
bottom layer, then, makes up the difference such that the total field capacity 
remains unchanged, giving a total soil depth that varies over the earth from 
50 to 225 cm. The lowest layer acts as a storage, supplying moisture to the 


83 



GISS Model II 


New 3-Laver Soil 

(After Abramopoulos et al., 1988) 


ET 

4 


4 

\ 

A 



n 



A 

No-Flux Boundary 

f 




Layers 1 and 2 
determined by 
field capacities 
specified by 
Matthews (1904) 


ET 


i 

it 

i i 

i 

4 

> 



A Layer 1 
1 always 10 cm 




A Layer 2 
X always 15 cm 

i L Depth of layer 
3 adjusted to 
maintain total 
field capacity of 
GISS Model II 

No-Flux Boundary 



/ / y~7^~7~7 / ~ 


- Transpiration only from 
top layer 

- Seasonal "on/off" toggle 
for exchange from lower 
layer to upper layer 

- Instantaneous moisture 
equilibration for vegetated 
areas 


- Solves diffusion equation 

- Transpiration from 

lower layers 


Figure 3.3. Comparison of new soil model to current GISS II 
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upper layers during dry periods. As shown by Entekhabi and Eagleson 

(1989b), the effectiveness of this mechanism is strongly dependent on soil 
texture (i.e., fractions of sand, silt, and clay). Heavier soils have a much 
greater ability to diffuse water upward against gravity due to high matric 
suction. Since this upward diffusion is replacing those features of GISS-II 

(i.e., prescribed growing season and instantaneous soil moisture equilibration 
between the two layers) designed to prevent summer desiccation, it will be of 
interest to examine summer continental temperatures, particularly where soil 
texture is composed mostly of sand. 

The two— layer heat conduction of the soil column remains the same as in 
GISS-II. However, the thickness of the top layer is made equal to the 
combined depths of the top two layers in the new soil moisture diffusion 
scheme. The bottom layers of both heat and moisture schemes are then of 
equal thickness. This matching provides a convenient common flux boundary 
over which heat transfer by movement of water can be computed along with 
heat conduction. Ice in the soil column contributes to relative soil saturation, 
but remains fixed in its soil layer of origin. 

In nature, the fraction of roots in different soil layers is dependent on 
many parameters including plant type, soil texture, and climatic region. It is 
not the purpose here to provide a full treatment of plant physiology in this 
parameterization. Rather, we sought to give a means of extracting water from 
lower soil layers to give a simple approximation to the function of plants in 
regional hydrologic balance. Based on observations compiled by Epstein 
(1973), it seems reasonable that for many plants, the root density in the soil 
decreases exponentially with depth. Seeking to find a reasonable value for root 
densities, the Screening Model was utilized to test sensitivity to various 
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weighting schemes. A weighting of 85% of plant roots in the top layer, 10% 
in the second, and 5% in the lowest layer was concluded to be a reasonable 
estimate. Higher percentages in the lower layers tended to dry the lowest 
layer in the screening model. We set these weightings uniformly 
over the globe since no firm data base exists for root density. 

In order to demonstrate the improvements brought about by the new soil 
diffusion scheme when compared with the current GISS scheme, we present the 
following plots of off-line sensitivity runs for three cases: 

1) A fifty-day response to an initially saturated top layer and initially 

dry bottom layer (s) with zero flux at the surface (Fig. 3.4a); 

2) A fifty-day day response to an initially saturated lower layer and 

initially dry top layer (s) with zero flux at the surface (Fig. 3.4b); 

3) A fifty-day response to a sinusoid forcing of the top layer (an 

implicit sinusoid flux at the surface) (Fig. 3.4c). 

A medium soil texture was used for the new soil diffusion scheme in 
these plots (the current GISS scheme has no dependence on soil texture). The 
thickness of layers follows the "nominal" thicknesses from Table 2.3 for the 
new soil diffusion scheme and an equivalent total thickness for the GISS 

diffusion scheme. Further, we are considering a vegetated area as opposed to 
bare soil (the vast majority of landsurface grids are predominantly vegetated). 

Figs. 3.4a-c show that, in general, the GISS scheme equilibrates at a 
much faster rate than the new scheme (this is true for other soil types as 
well). In fact, in the case of upward diffusion (Fig. 3.4b) instantaneous 
equilibration of relative soil saturation in the GISS scheme is evident. The 
equilibrium state for the GISS scheme, regardless of the relative thicknesses of 
the two layers, is for the relative soil saturation values of all layers to be 
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Figure 3.4a 


GISS Diffusion: 

Top Layer Initially Saturated 



New Soil Diffusion: 

Top Layer Initially Saturated 



Fifty day off-line comparison of GISS (top) and New Soil 
Diffusion (bottom) for an initially saturated top layer and dry 
lower layer (New Soil Diffusion uses medium soil texture, 
GISS scheme has no dependence on soil texture) 
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GISS Diffusion: 

Bottom Layer Initially Saturated 



New Soil Diffusion: 
Bottom Layer Initially Saturated 



Figure 3.4b 


Fifty day off-line comparison of GISS (top) and New Soil 
Diffusion (bottom) for an initially saturated lower layer and dry 
upper layer (New Soil Diffusion uses medium soil texture, GISS 
scheme has no dependence on soil texture) 
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GISS Diffusion: 

Sinusoidal Forcing of Top Layer 



New Soil Diffusion: 
Sinusoidal Forcing of Top Layer 



Day 


Figure 3.4c Fifty day off-line comparison of GISS (top) and New Soil 
Diffusion (bottom) for a sinusoidal forcing of upper layer (New 
Soil Diffusionuses medium soil texture; GISS scheme has no 
dependence on soil texture) 
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equal. This is unrealistic since gravity ought to be forcing a higher 
concentration in the lower layer(s) in the absence of other forcings. In the 
new soil scheme, the effects of gravity are accounted for and we see in 
Figs. 3.4a-c that indeed the lower layers have preferentially higher equilibrium 
relative saturations. Finally, in the sinusoidal forcing case particularly, we see 
that the response of the lowest layer to the upper layer is very dampened in 
the new scheme as compared to the GISS scheme. The quick response in the 
GISS scheme is at least to some degree a surrogate for plant transpiration. In 
the new scheme, transpiration by root extraction is modeled directly, although 
these six plots have not included these effects since it would require 
complications which would cloud the messages of the plots. The rate at which 
diffusion takes place between the layers will have significant influence on the 
model climate. The amplitude of the annual cycle in surface soil saturation 
will affect the amplitude of surface heat and moisture balances among other 
factors. These effects will become evident when the new Abramopoulos et al. 
(1988) scheme for soil moisture diffusion is implemented in the GISS-GCM 
model and tested along with the new landsurface hydrology parameterization 
with subgrid scale spatial variability. 

C. Results 

We have discussed so far the changes to be explored within the 
landsurface hydrological parameterization of the GISS GCM. From here on we 
discuss the results as obtained from the simulations with the three-dimensional 
GISS GCM II. 
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1 Fundamental Changes in Hydrology Over Landsurface Grid Squares 

The changes induced by the new hydrology parameterization may first be 
seen by examination of individual landsurface grids. We consider here diag- 
nostics of single points in space averaged over increments of time as well as 
distributions of diagnostics of all landsurface grids averaged over the entire 
simulation period. More attention is given to the results of Run G— 3 (see 
Table 3.1) since it produced the strongest changes. 

The runoff ratio (the ratio of runoff to precipitation) at a point may be 
averaged over time and plotted against the average relative surface soil 
moisture. In this way we see the actual model results which compare to the 
off-line analytical plots presented earlier. Fig. 3.5 shows monthly-averaged 

runoff ratios in both the Control Run G-0, and Run G-3 over a Central 
Argentina grid square. This particular grid experienced a large seasonal range 
of soil moisture values. It is seen that the plot for the Control Run G— 0 
deviates slightly from the empirical R = js line which is inherent in the 
model. This is due to ponding on the surface (giving a slightly higher value 
of R) and also to the surface storage of snow (giving a slightly lower value of 
R). The GISS scheme has no dependence on precipitation intensity, soil type, 
or spatial characteristics, but only on relative surface soil moisture. 

The bottom plot of Fig. 3.5 shows the runoff ratio of the same grid 
square, but for Run G-3. Immediately obvious is the nonlinear relationship of 
R with s. High values of s result in runoff ratios reaching 0.7 in the monthly 
mean. The current GISS runoff parameterization would not allow such high 
values under any conditions except a persistently saturated surface layer. The 
effects of the higher runoff ratios in general have caused a higher likelihood of 
lower values of relative soil saturation (evident near the left of the graph). 
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Figure 3.5 Comparison of Control Run G-0 and Run G-3 (space/ 
soil /storm) of the monthly runoff ratio over a Central Argentina 
gria (Note: Deviations from R = 1/2 s line in Control Run are 

due to ponding on the surface and snow storage) 
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Actual soil moisture data is unavailable for this region (as it is in general) 
and comparisons against "ground truth" are not possible. The important point 
is that the runoff ratio in Run G-3, with the effects of spatial variability and 
realistic infiltration dynamics, exhibits a nonlinear relationship with soil 
saturation. 

When various hydrologic fluxes for all landsurface grids are plotted, 
patterns of direction and magnitude of changes may be seen across the globe. 
Figs. 3.6 through 3.8 illustrate the difference between the G-3 Run and 
Control Run of the diagnostics of runoff and evaporation, plotted against 
relative soil moisture and precipitation intensity. Figure 3.9 gives the 
distribution of soil saturation versus precipitation intensity. 

Figs. 3.6 and 3.9 show that grids with low soil saturation are often 
characterized by high precipitation intensity as well as high runoff ratio. Since 
the offline plots of R (Figs. 3.1a-c) show greater sensitivity at high rainfall 
rates, it is expected that the largest difference between Runs G-0 and G-3 are 
at the low soil saturations (see Figs. 3.6 and 3.8): the strongest positive 

changes occurred in grids having low values of s. 

When plotted against precipitation intensity (Fig. 3.7) the runoff 
difference shows a very clear envelope rising sharply at low intensities and 
leveling off as intensities exceed 25 mm/day. Based on the off-line analysis 
cliscussed earlier, one might have expected the slope of this envelope to 
continue to be quite steep even as the precipitation intensities became high. 
However, in high intensity rainfall areas, the current GISS model produces 
generally high soil moistures, which give rise to higher runoff. In Run G-3, 
these soil moistures have generally decreased, and therefore the effects of 
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Figure 3.6 


Grid-by-grid plot of the difference in runoff produced in 
Run G-3 (space/soil/storm) and the Control Run G-0, plotted 
against relative surface soil moisture for both cases 
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Figure 3.7 Grid-by-grid plot of the difference in runoff produced in 
Run G-3 (space/soil/storm) and the Control Run G-0, plotted 
against precipitation intensity for both cases 
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Figure 3.8. 
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Grid-by-grid plot of the difference in actual evaporation 
produced in Run G-3 (space/soil/storm) and the Control Run 
G-0, plotted against relative soil saturation for both cases 
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Surface Relative Soil Moisture vs. 
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Figure 3.9 Distribution of relative surface soil saturation and precipitation 
intensity for the Control Run G-0 and Run G— 3 

(space/ soil /storm) 






spatial variability, though still strong, are held in check by this negative 
feedback mechanism operating in high intensity rainfall regimes. 

In the case of evaporation, we see in Fig. 3.8 that the strongest effects 
are also in grids with low relative soil saturation. Evaporation changes in 
Run G-3 were almost everywhere negative. This would appear to be due to 
the drying of soils produced by increased surface runoff. Because potential 
evaporation is in general extremely high (a point we shall return to later) 
evaporation is very efficient in both the GISS scheme and the spatial 

variability scheme. The main reason for the negative differences at low values 
of s in Fig. 3.8 is that there is simply not as much water available to 
evaporate since a greater portion has been lost to runoff. The underestimation 
of continental runoff in much of the 8' * 10° GISS GCM model is coupled 
with the severe overestimation of land surface evaporation in that model. 
Simulation Run G-3, which includes the effects of subgrid scale spatial 
variability on the estimation of surface hydrologic fluxes, significantly improves 
the model climate runoff rates over the continents; as a consequence, the 
remainder of the water balance and components such as the land surface 

evaporation rates are also improved with respect to correspondence to 
observations. 

Next we consider the effects of the new soil moisture diffusion scheme. 
Figs. 3.10a-b and 3.11a-b show the time series of relative soil moistures of a 
Northwest U.S. grid (mainly covering Montana, chosen for the large range of 
soil moisture values) and an eastern Sahel grid (a dry desert area). The top 
soil layer (layer 1) relative soil saturation has been collected and plotted as a 
special daily diagnostic whereas the lower layer is the standard monthly- 

average diagnostic. Note that for Runs G-2 and G-3, which actually used 
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[Control Run (G-0)] 



Northwest U.S. Grid Relative Soil Saturation 
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Figure 3.10a Three-year plots of daily averaged soil saturation in layer 1 and 
monthly averaged soil saturation in layer 2 for Control Run G-0 
and Run G-l (space) over a Northwest U.S. grid 
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Northwest U.S. Grid Relative Soil Saturation 
[Run G-2 (sp/so)] 



Northwest U.S. Grid Relative Soil Saturation 
[Run G-3 (sp/so/st)] 



Figure 3.10b Three-year plots of daily averaged soil saturation in layer 1 and 
monthly averaged soil saturation in Runs G-2 (space/soil) and 
G-3 (space/soil/ storm) over a Northwest U.S. grid 
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Figure 3.11a Three-year plots of daily averaged soil saturation in layer 1 and 
monthly averaged soil saturation in layer 2 for Control Run G-0 
and Run G-l (space) over an Eastern Sahel grid 
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Eastern Sahel Grid Relative Soil Saturation 
[Run G-2 (sp/so)] 



Eastern Sahel Grid Relative Soil Saturation 
[Run G-3 (sp/so/st)] 



Figure 3.11b 


Three-year plots of daily averaged soil saturation in layer 
monthly averaged soil saturation in layer 2 for Run 
(space/soil) and G-3 (space/soil/storm) over an Eastern 

grid 


1 and 
G-2 
Sahel 
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three soil layers, the "layer 2" designation represents both layers 2 and 3. 
This is done for ease of comparison since the field capacity of layer 2 in the 
GISS scheme is roughly equivalent to the combined field capacities of layers 2 
and 3 in the new 3-layer scheme. 

As expected, there is only a slight change between the Control Run and 
Run G-l (which retained the same GISS soil moisture diffusion) in the U.S. 
(Fig. 3.10a). The seasonal cycle has generally the same shape, and the 
seasonal range is the complete range from 0 to 1. It should be noted that the 
G-l Run was started on July 1 instead of January 1, creating a half-year lag 
in its cycle. The G-l Run is slightly more moist, due to a decrease in runoff 
in this particular grid. In both of these plots (Fig. 3.10a) we see the effects 
of the GISS instantaneous upward diffusion for the drying phase of layer 2 
over vegetated areas. This is evident by the closeness of the drying curves of 
layers 1 and 2 (note the far greater independence of these curves for the 
wetting phases). This grid is not completely vegetated. If it were, the curves 
for layers 1 and 2 would match identically during the drying phase of layer 2 
(as seen previously in the off-line diffusion Figs. 3.4b and c). Although not 
obvious from these plots, there is also a prescribed growing season for 
vegetated areas running from May 1 to September 1 in the Northern 

Hemisphere in mid- to high latitudes. Outside this time period in the GISS 
scheme, layer 2 soil moisture can only be increased, not depleted. Then it is 
allowed to supply water to layer 1 during the growing season. 

For the case of Run G-2, the new 3-layer soil column is seen to yield 

very similar results for layer 1: the shape and magnitude of the seasonal 

cycle are roughly the same (Fig. 3.10b). While the amplitude of layer 2 soil 
moisture is also roughly the same, its mean value is about 30% higher. 
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Apparently the physically based equations of soil water movement are able to 
supply water efficiently to the point that the artificial prescribed growing 
season and unrealistic instantaneous upward diffusion over vegetated areas are 
no longer necessary. There is a much greater independence of the two layers, 
and there is now a consistent 1- to 2-month lag between the two cycles, 
whereas before there was a 1— to 2— month lag only at the peak of soil 
moisture (in winter) and almost no lag in the drying phase and at the low 
point at the end of summer. 

Finally, in Run G-3 with the increased surface runoff, all three layers 
are consistently drier. The most significant change over the U.S. grid is the 

prolonged dry period in the summer months, which appears to be 

approximately twice as long as the dry period in the other runs. 

Next we consider the eastern Sahel grid in Figs. 3.11a-b. Since this is a 
completely non— vegetated, low latitude grid, there are no instantaneous upward 
diffusion effects or growing season effects operating in the GISS scheme. Thus 
we see a lag in both the wetting and drying phases of the Control Run (G-€). 
Rain is relatively infrequent and therefore it is possible to see the intermittent 
pattern of completely dry soil to wet, and back to dry again from one rainfall 
event. In Run G-l (Fig. 3.11a), we see first of all that the soil is more 

moist than the Control Run. This is a direct result of a lower runoff 
coefficient for light soil and low precipitation intensity (see upper plots of 
Figs. 3.1a and 3.1b) and also a lower bare soil evaporation efficiency. 

Reductions in runoff and evaporation allow for more moisture in the soil. 

When the new soil diffusion scheme is incorporated (Run G-2, top of 
Fig. 3.11b) there is a fundamental difference: the lower layer response to 
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changes in the upper layer is severely damped. In fact, even over long periods 
of time without rain the lower layer is relatively unaffected. The new runoff 
parameterization allows more infiltration of greater intensity rainfall events, 
thus the spikes of upper layer soil moisture are only half as high on average 
as those of the GISS soil diffusion. Interestingly, the time constants of the 
wetting and drying phases of the upper layer appear to be unchanged. 

2. Global Hydrologic Balance 

a. Globally Averaged Water Balance . We now turn to the analysis of 
the simulation results compared with observed data. We follow a progression 
from the very largest scale comparisons (globally averaged quantities) to more 
detailed (zonally averaged and continentally averaged quantities) and finally to 
mapped distributions over the earth. The mapped distributions are found in 
Appendix A. 

The most basic picture of hydrologic performance in the model is the 
globally averaged hydrologic cycle. There are several sources of compiled 
values for the globally averaged precipitation and evaporation over the oceans 
and landsurface, and the exchange of moisture between those two regions 
(Budyko, 1963; Baumgartner and Reichel, 1975; Budyko, 1978; Henning, 1989; 
and others). In the long-term mean, of course, this latter exchange must 
equal the runoff from the land. As a comparison, we have plotted schematics 
of the global hydrologic balance of all simulations along with one set of 
observed values (Budyko, 1978) in Fig. 3.12. The values are all in annual 
depths [centimeters] per unit area, and the two exchange values (atmospheric 
divergence over land and runoff) are normalized by total landsurface area. 
Run G-3 is slightly out of balance due to incomplete drainage of initial 
storage from the landsurface. Because of this, it is expected that the runoff 
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Figure 3.12. Global hydrologic balance [in cm per year] for all simulations 
with observations from Budyko (1978) 


106 









value in Run G-3 is slightly higher than what would result from a longer 
duration simulation. The ocean evaporation value is essentially constant for 
all runs due to the use of prescribed ocean surface temperatures. 

First, we notice the GISS II model (Control Run) has underestimated 
runoff and overestimated both precipitation and evaporation over the 
landsurface (Fig. 3.12). The range of values of global landsurface runoff found 
by the observers cited earlier is from 27 to 35 cm annually. The GISS II 
value does not even fall in this range. Further, neither the GISS II 
precipitation over land nor evaporation over land fall into the ranges of the 
observers for these quantities (precipitation over land observations range from 
72 to 80 cm annually whereas GISS II gives 88 cm; evaporation over land 
observations range from 41 to 48 cm annually whereas GISS II generates 65, 
roughly 50% more than observed). Similarly, Runs G— 1 and G-2 grossly 
overestimate landsurface evaporation and underestimate runoff. 

Major improvements occurred with Run G-3 (with reduced fractional 
wetting parameter of moist convective events; from 0.6 as in Run G-2 to 
0.15). Here we see precipitation over land and evaporation over land falling 
to near their observed ranges. Global runoff is slightly high, but as mentioned 
earlier the runoff value will likely decrease as soil storages are depleted to 
their equilibrium states. Agreement with the Budyko (1978) values is 

remarkable. It is interesting to note the effects of the landsurface hydrology 
on the precipitation over oceans. The lower ocean precipitation value, while 
lower than Budyko (1978), is still well within the range of the other observers 
(105-127 cm annually). 

Inside the schematic landsurface boxes in Fig. 3.17 is the average 
landsurface temperature. Observed values are derived from Legates and 
Willmott (1990). In Hansen et al. (1983) it was noted that summertime 
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temperatures over continents were 2—3° C too cool in GISS II. As can be 
seen in the schematic, landsurface temperatures rose 0.7° C on average in 
Run G— 3. As will be seen on the global temperature maps, this warming 
occurred mainly in the tropics over the Amazon, Central America and the 
Caribbean as well as North Africa in June-July-August, and over northern 
Australia in December-January-February. Aside from this slight warming, 
though, the G-3 Run has shown a significant improvement in the 
globally-averaged hydrology of the GISS model. 

h Tonally Averaged Water Balance, Figs. 3.13a-c show the zonally 
averaged values of precipitation, evaporation, and surface runoff over the 
landsurface for each of the simulations as well as observed values compiled by 
Zubenok (1970). Although Hansen et al. (1983) plotted zonally-averaged 
precipitation and evaporation and found very good agreement with 
observations, their comparison included the oceans. Since two-thirds of the 
earth is covered by oceans and the hydrologic fluxes over oceans are 
considerably larger than over land, the contribution of the landsurface alone 
was masked. Here we see that the zonally averaged hydrology over 
landsurfaces in the GISS Model II is actually in rather poor agreement with 
data (Run G-0). It should be noted that in zonally averaged plots, each 
latitude belt is given equal weight. In reality, however, the actual volume of 
the fluxes depends on the amount of landsurface in the latitude belt. Since 
most of the land mass over the earth is in the Northern Hemisphere (see 
Fig. 3.14), we must keep in mind that errors in zonal plots should be 
interpreted accordingly. In latitude bands where land comprises a small 
fraction of the total area, the estimates will have larger variance. 
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Figure 3.13a Zonally averaged annual landsurface precipitation 
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Figure 3.13b Zonally averaged 
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Figure 3.13c Zonally averaged annual surface runoff 
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latitude (taken from the GISS GCM) 


function of 


112 


As presented here, we see that in general, precipitation (Fig. 3.13a) is 
modeled quite well in the tropics, but agreement with observations deteriorates 
in the mid- and high latitudes. Particularly in the Northern Hemisphere 
mid-latitudes (where data are generally more readily available) precipitation is 
overestimated. The only strong deviation from the Control Run is Run G-3 
in the tropics. Here we see that although Run G-3 had a very good globally 
averaged precipitation over land, it is actually underestimating in the tropics 
and overestimating in mid- and high latitudes. 

In Fig. 3.13b the simulations G-l and G-2 are essentially equivalent to 
the control in the zonal evaporation over land. Run G-3, however, gives 
agreement with observations in the tropics which is remarkable. In Fig. 3.13c 
the same pattern holds: runoff in Run G-3 comes much closer to the 

observed values in the tropics. The spike at 50° south latitude in the 

observed runoff is due to the fact that there is very little land in this belt 
(Fig. 3.14). Errors in data-sparse regimes such as this may be reflecting 
measurement error. Since evaporation observations are generally the residual 
of precipitation and runoff, an error in runoff measurement necessitates an 
error in evaporation "observations". 

The balance at the tropics in Run G-3 makes it clear that the model is 
simply not producing enough precipitation in this region. Obviously if the 
evaporation is virtually in agreement with observation and both precipitation 
and runoff are low, there is a problem not in the landsurface parameterization, 
but in the atmospheric moisture transfer to the tropics and/or in the moist 
convection parameterization. 

Conversely, in the Northern Hemisphere mid— latitudes, where all three of 
these variables are too high, there is an excess of precipitation. This again 
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cannot be solved by the landsurface parameterization since runoff and 
evaporation are also exceeding observations. 

We have seen that the effects of the reduced fractional wetting in moist 
convective rainfall in Run G-3 are strongest by far in the tropics. This 
warm, high humidity region has both the highest amount of rainfall and the 
highest percentage of moist— convective rainfall (see Fig. 3.2). Moving away 
from the tropics, the effects in Run G-3 are reduced nonlinearly, preserving 
similar values to the GISS Model II. As seen in these zonally averaged 
landsurface figures, the results are quite good for the GISS GCM equipped 
with the Entekhabi and Eagleson (1989a) new land surface hydrology 
parameterization with subgrid scale spatial variability and the new 
Abramopoulos et al. (1988) soil diffusion scheme. We wish to reemphasize in 
particular the dramatic improvement of the evaporation on all scales global, 
zonal, and, as will be seen, in continental as well as over the river basins. 

r. Continental Water Balance. As the next step in our progression to 
finer spatial detail, we consider the water balance over the continents. Shown 
in Table 3.3 and Figs. 3.15a-c are the precipitation, evaporation and surface 
runoff over continents (excluding Antarctica) for the four simulations as well 
as observations by Henning (1989). North America does not include 
Greenland, and Australia does not include New Zealand. 

As we have seen before, Runs G-l and G-2 differ only slightly from the 
control in most cases. Some modest exceptions include the substantial increase 
in both precipitation and evaporation over Europe in Run G-l (which changed 
runoff very little), and a significant decrease in runoff over Africa in both G-l 
and G-2. In Run G-3, there are major changes on all fronts. On all 
continents runoff was increased, while precipitation and evaporation were 
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Continental Hydrologic Balance (m/year) 


Continent 

South America 
Africa 
Australia 
Asia 

North America 
Europe 


Continent 

South America 
Africa 
Australia 
Asia 

North America 
Europe 


Continent 

South America 
Africa 
Australia 
Asia 

North America 
Europe 


i. VV1J 

G-0 

^ -L VJ. 

G-l 

VII till 

G-2 

11IV Y CCU i 

G-3 

Observations 

(control) (sp) 

(sp/so) 

(sp/so/st) (Henniner. 1989) 

1396 

1432 

1361 

1178 

1481 

1017 

1003 

911 

831 

695 

768 

734 

756 

597 

471 

1016 

1064 

1052 

896 

699 

778 

851 

796 

734 

622 

767 

911 

865 

717 

597 

Evanoration (mm/vear) 


G-0 

G-l 

G-2 

G-3 

Observations 

(control) 

[££} 

(sd/sq) (sn/so/st) (Hennine. 1989) 

1172 

1242 

1172 

682 

883 

877 

902 

828 

478 

502 

662 

672 

674 

400 

403 

708 

776 

737 

519 

417 

508 

594 

549 

437 

394 

562 

695 

646 

491 

362 

Surface Runoff (mm/vear) 


G-0 

G-l 

G-2 

G-3 

Observations 

(control) 

(sp) 

(sp/so) (sp/so/st) (Henning. 1 989) 

216 

180 

185 

516 

598 

139 

100 

89 

373 

193 

105 

65 

82 

208 

67 

309 

291 

309 

409 

282 

270 

255 

244 

326 

228 

198 

214 

208 

243 

235 


Table 3.3. Hydrologic balance of continents with observations by Henning 
(1989) 
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Figure 3.15a. Annual continental precipitation 
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Evaporation (mm/year) 
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Figure 3.15b. Annual continental evaporation 
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Runoff (mm/year) 


Continental Runoff 
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Figure 3.15c. Annual continental runoff 
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decreased. By comparison to observations, this was a significant improvement 
for all continents for evaporation, and an improvement in all but 
South America for precipitation. However, the runoff was only improved for 
two continents (South America and Africa) while the other four worsened in 
the direction of overestimation. It should be noted that the estimates of 
evaporation and runoff in Henning (1989) are independent of that made by 
Zubenok (1970). According to both sources, the Control Run G-0 with the 
GISS GCM severely overestimates landsurface evaporation and underestimates 
runoff. Hydrologic fluxes in Run G-3 with the improved parameterization are, 
however, in agreement with both sources. 

The problem with runoff over continents is not solely a result of the 
landsurface parameterization but is partly due to there being an excess of 
precipitation. This can be seen by first noting that Run G-3 evaporation 
values over continents with high runoff values are all either very close to 
observed or else too high. But in all cases, the overestimation of precipitation 
is still higher, and thus, since runoff is also too high, we see that as before 
there is simply too much water being fed into the system by the atmospheric 
transport and rainfall generation. If the evaporative input to the atmosphere 
is roughly correct, then the problem of excessive precipitation lies in the 
atmospheric branch of the hydrologic cycle, which is not the subject of this 
work. 

d. Major River Basin Water Balances. As a further step in evaluating 
the model hydrologic performance, the runoff over major river drainage systems 
will be compared to observations. We follow Russell and Miller (1990) in 
discretizing the GISS model into 33 river basins around the earth. In their 
paper, a 4° x 5° grid spacing was used in the model, as opposed to an 
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8* * 10’ spacing used here. Further, their results were plotted by volume, 
whereas we will plot per unit area values in the water balance. The reason 
for this is that the modifications made act on a per unit area basis over grid 
squares and the pattern of the effects of the modifications is more readily 
apparent in this form. Furthermore ranking according to per unit area runoff 
helps in isolating the changes according to the levels of humidity in climate. 

Tables 3.4a-c give the water balance components of the 33 basins in a 
diverse range of climatic zones arranged in descending order of per unit area 
observations. Figs. 3.16a-f show the runoff plotted in descending order of per 
unit area runoff; the first river basins shown are thus generally the most 
humid ones and the ones with the highest runoff ratios. It can be seen again 
that precipitation in the Northern Hemisphere is rather consistently 
overestimated in most of the basins. 

In some cases, such as the Yellow River, the overestimation of 
precipitation is severe. The effects of this may be seen in the gross 

overestimation of runoff for this river (more than 15 times the observed 
value), in spite of evaporation being twice as large as observations. The 
observed value of runoff for this narrow basin is also prone to underestimation 
due to overbank flooding at times of high flow. 

The striking improvement due to Run G-3 apparent in the tropical river 
basins depicted in Figure 3.16a-f confirms what was seen earlier in the 
zonally— averaged runoff plot for tropical regions. As the total runoff value 
decreases in more arid basins, we see the agreement becoming worse. 
However, one must keep in mind the overestimation of precipitation in these 
basins. 
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Major River Basin Precipitation (m/year) 


River 

G-0 

G-l 

G-2 

G-3 

Observations 

Basin 

(control) (sp) 

(sd/so) (sD/so/st) (Russell/Miller. 1990) 

Irrawady 

2.07 

2.30 

2.19 

1.90 

1.92 

Amazon 

1.53 

1.50 

1.48 

1.43 

1.90 

Magdalena 

2.17 

2.12 

2.14 

1.81 

1.55 

Orinoco 

2.35 

2.29 

2.31 

2.08 

1.55 

Congo 

1.84 

1.74 

1.55 

1.38 

1.47 

Hsi Chiang 

1.73 

1.66 

2.11 

1.41 

1.43 

Mekong 

2.27 

2.43 

2.41 

2.01 

1.38 

San Francisco 

2.29 

2.38 

2.07 

1.66 

1.36 

Brahma-Ganges 

1.05 

1.41 

1.38 

0.99 

1.16 

LaPlata 

0.62 

0.64 

0.65 

0.60 

1.15 

Zambesi 

1.30 

1.18 

1.15 

1.12 

1.01 

Yangtze 

1.77 

1.73 

1.84 

1.56 

1.00 

Niger 

1.08 

1.16 

1.00 

0.70 

1.00 

St. Lawrence 

1.05 

1.09 

0.96 

0.89 

0.88 

Fraser 

0.99 

1.20 

1.04 

0.94 

0.86 

Mississippi 

0.90 

0.93 

0.92 

0.84 

0.75 

Danube 

0.92 

0.99 

0.94 

0.83 

0.70 

Nile 

1.23 

1.12 

1.00 

1.13 

0.68 

Columbia 

0.88 

1.13 

0.94 

0.86 

0.63 

Indus 

0.82 

0.78 

0.66 

0.65 

0.55 

Murray 

0.92 

0.88 

0.93 

0.82 

0.54 

Sevemay Dvina 

0.60 

0.75 

0.73 

0.60 

0.51 

Amur 

0.96 

1.01 

0.99 

0.94 

0.50 

Yellow 

2.10 

1.97 

2.04 

2.07 

0.48 

Yukon 

0.84 

0.91 

0.89 

0.90 

0.44 

Ob 

0.34 

0.41 

0.36 

0.33 

0.42 

Orange 

0.64 

0.63 

0.60 

0.65 

0.40 

Tigris-Euphrates 

0.49 

0.61 

0.59 

0.53 

0.39 

McKenzie 

0.68 

0.74 

0.71 

0.67 

0.39 

Yenesei 

0.53 

0.57 

0.54 

0.47 

0.37 

Lena 

0.57 

0.65 

0.62 

0.58 

0.33 

Kolyma 

0.72 

0.74 

0.74 

0.72 

0.26 

Colorado 

0.85 

1.28 

1.02 

1.06 

0.25 


Table 3.4a Per unit area annual precipitation over major river basins; 

observations compiled by Miller and Russell (1990) (boldfaced 
lines indicate to five rivers by volume ) 
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Major River Basin Evaporation ( m / year ) 


River 

G -0 

G-l 

G -2 

G -3 

Observations 

Basin 

(control) 

(sp) 

(sp/so) (sp/so/st) 

(Russell/Miller. 1990) 

San Francisco 

1.76 

1.82 

1.65 

0.97 

1.21 

Congo 

1.54 

1.58 

1.43 

0.77 

1.14 

LaPlata 

0.57 

0.61 

0.60 

0.38 

0.99 

Irrawady 

1.62 

1.84 

1.74 

1.19 

0.93 

Amazon 

1.33 

1.33 

1.28 

0.76 

0.88 

Niger 

0.98 

1.08 

0.92 

0.36 

0.85 

Zambesi 

1.12 

1.09 

1.05 

0.66 

0.83 

Mekong 

1.68 

1.88 

1.82 

1.21 

0.78 

Hsi Chiang 

1.28 

1.34 

1.49 

0.88 

0.74 

Nile 

1.05 

1.00 

0.91 

0.61 

0.65 

Mississippi 

0.67 

0.74 

0.74 

0.60 

0.57 

Magdalena 

1.79 

1.89 

1.81 

0.95 

0.56 

Yangtze 

1.11 

1.20 

1.15 

0.86 

0.54 

Murray 

0.77 

0.80 

0.83 

0.58 

0.52 

Brahma-Ganges 

0.86 

1.15 

1.10 

0.70 

0.50 

St. Lawrence 

0.62 

0.75 

0.64 

0.55 

0.45 

Orinoco 

1.91 

1.99 

1.94 

1.08 

0.44 

Danube 

0.71 

0.82 

0.77 

0.62 

0.44 

Yellow 

0.89 

0.86 

0.83 

0.81 

0.42 

Orange 

0.57 

0.57 

0.56 

0.46 

0.39 

Tigris-Euphrates 

0.44 

0.56 

0.56 

0.41 

0.35 

Fraser 

0.64 

0.85 

0.73 

0.66 

0.35 

Amur 

0.64 

0.69 

0.65 

0.58 

0.32 

Indus 

0.57 

0.61 

0.54 

0.42 

0.30 

Ob 

0.24 

0.32 

0.26 

0.22 

0.27 

Columbia 

0.67 

0.88 

0.76 

0.67 

0.26 

Colorado 

0.65 

1.02 

0.86 

0.77 

0.22 

McKenzie 

0.39 

0.45 

0.44 

0.39 

0.22 

Sevemay Dvina 

0.31 

0.41 

0.38 

0.30 

0.21 

Yukon 

0.36 

0.43 

0.40 

0.37 

0.20 

Kolyma 

0.28 

0.29 

0.31 

0.26 

0.15 

Yenesei 

0.38 

0.44 

0.39 

0.33 

0.15 

Lena 

0.39 

0.43 

0.42 

0.37 

0.13 


Table 3.4b Per unit area annual evaporation over major river basins; 

observations compiled by Miller and Russell (1990) (boldfaced 
lines indicate top five rivers by volume ) 


122 



Major River Basin Runoff (m/year) 


River 

Basin 


G-0 G-l 
(control) (sp) 


G-2 G-3 Observations 



Orinoco 

0.43 

0.29 

0.34 

1.03 

1.11 

Amazon 

0.20 

0.17 

0.20 

0.70 

1.02 

Irrawady 

0.41 

0.45 

0.43 

0.83 

1.00 

Magdalena 

0.36 

0.24 

0.30 

0.89 

0.99 

Hsi Chiang 

0.44 

0.33 

0.58 

0.64 

0.69 

Brahma-Ganges 

0.20 

0.27 

0.25 

0.31 

0.66 

Mekong 

0.55 

0.54 

0.57 

0.92 

0.59 

Fraser 

0.35 

0.36 

0.32 

0.30 

0.51 

Yangtze 

0.63 

0.54 

0.68 

0.78 

0.46 

St. Lawrence 

0.43 

0.35 

0.31 

0.37 

0.43 

Columbia 

0.20 

0.25 

0.21 

0.24 

0.37 

Congo 

0.29 

0.16 

0.13 

0.65 

0.33 

Sevemay Dvina 

0.28 

0.33 

0.34 

0.32 

0.30 

Indus 

0.25 

0.17 

0.10 

0.25 

0.25 

Danube 

0.20 

0.18 

0.17 

0.24 

0.25 

Yukon 

0.47 

0.48 

0.48 

0.52 

0.23 

Yenesei 

0.14 

0.14 

0.14 

0.17 

0.22 

Lena 

0.18 

0.21 

0.20 

0.24 

0.21 

Zambesi 

0.19 

0.09 

0.12 

0.49 

0.19 

Amur 

0.34 

0.33 

0.32 

0.38 

0.18 

Mississippi 

0.22 

0.19 

0.18 

0.27 

0.18 

LaPlata 

0.05 

0.04 

0.05 

0.23 

0.17 

McKenzie 

0.29 

0.28 

0.27 

0.29 

0.17 

Niger 

0.10 

0.09 

0.08 

0.35 

0.16 

San Francisco 

0.53 

0.55 

0.43 

0.75 

0.15 

Ob 

0.10 

0.10 

0.11 

0.12 

0.15 

Kolyma 

0.41 

0.45 

0.43 

0.46 

0.11 

Yellow 

1.21 

1.11 

1.20 

1.26 

0.07 

Tigris-Euphrates 

0.06 

0.05 

0.03 

0.12 

0.04 

Colorado 

0.19 

0.27 

0.14 

0.28 

0.03 

Nile 

0.17 

0.12 

0.10 

0.54 

0.03 

Murray 

0.15 

0.08 

0.11 

0.25 

0.02 

Orange 

0.07 

0.05 

0.04 

0.20 

0.01 


Table 3.4c. Per unit area annual runoff over major river basins; observations 
compiled by Miller and Russell (1990) (boldfaced lines indicate 
top five rivers by volume) 
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Figure 3.16a Per unit area annual runoff over major river basins 
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Figure 3.16b Per unit area annual runoff over major river basins 
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Figure 3.16c Per unit area annual runoff over major river basins 
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Figure 3.16d Per unit area annual runoff over major river basins 
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Figure 3.16e Per unit area annual runoff over major river basins 
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Figure 3.16f Per unit area annual runoff over major river basins 
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3 Temperature and the Heat Balance 

The hydrologic balance over landsurfaces has a strong effect on the 
regional heat balance. The mechanisms include: 1) atmospheric transport of 

latent heat, 2) the radiative absorption and emittance by atmospheric water 
vapor, 3) the reflective properties of clouds, 4) the reflective properties of 
snow on the earth’s surface, 5) the increase of soil heat capacity due to the 
presence of soil moisture and 6) the release of latent heat from the surface by 
evaporation. There are also secondary effects such as those due to the 
dependence of vegetation type and productivity on the moisture budget and 
the influence of the heat budget on wind patterns. Of all these mechanisms, 
it is the evaporation which has the strongest effect on the surface temperature; 
it is a major cooling mechanism of the landsurface when moisture is available. 
As moisture in the soil is depleted, the heat partitioning at the surface shifts 
toward higher rates of sensible heat flux. Since the sensible heat flux is a 
much less efficient transfer process, temperatures at the surface may rise 

significantly under reductions in soil moisture. 

In Appendix A, Figs. A.5a-d show the global distribution of surface air 
temperature of the Control Run G-0 for December-January-February (DJF) 
and June-July-August (JJA), as well as deviations from the Control for each 
simulation G-l to G-3. While there are isolated changes (mainly summer 
cooling) in Runs G-l and G-2, the most dramatic effects are in Run G-3 
since the hydrology changed significantly in this simulation. By comparison to 
the evaporation Figs. (A.2a-d) and surface soil moisture Figs. (A.4a-d) it can 
be seen that generally lower temperatures are produced in areas of increased 
soil moisture and increased evaporation. The converse is also true. Our main 
concern is with Run G-3 since it showed the greatest improvement 
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hydrologically. We now consider surface air temperature comparisons with 
data in zonal mean and global distribution. This parameter captures the 
effects of the new schemes on the heat balance in the GCM most effectively. 

In the zonal mean, the effects of the modified hydrology become more 
apparent. Figs. 3.17a-c show zonal mean surface air temperatures over the 
landsurface for all simulations along with observations compiled by Legates and 
Willmott (1990). Fig. 3.17a is the annual mean while Figs. 3.17b and 3.17c 
are the winter and summer means, respectively. As can be seen, all 
simulations are essentially the same except in the tropics. Here Run G-3 
diverges from the other three, ranging from 0 to 2 degrees warmer. In the 
annual mean, we see that Run G— 3 is actually in exceptionally good agreement 
with data from 20 S. Latitude to 20 N. Latitude. The seasonal means also 
show Run G-3 to be an improvement except for a small latitude belt near the 
equator in the summertime. In the global distribution, we shall see this is 
mainly caused by excessive heating in the Amazon and Caribbean regions. 

Fig. 3.18 gives observed January and July global surface air temperature 
distributions. Figs. 3.19a and b show the Control Run G-0 and Run G-3 
DJF and JJA surface air temperature distributions. Due to reduced runoff, 
and therefore reduced soil moisture and evaporation, Run G-3 is in general 
slightly warmer than the Control Run. For the winter months, agreement is 
quite good. For the summer, Run G-3 has the Caribbean and northern 
South America to be approximately 2-4 degrees too warm. Also, the region 
over Saudi Arabia has gotten slightly warmer (1 degree), but as it was 
already on the order of 3-5 degrees too warm in the Control Run, this is not 
severe. The region over west Africa has improved, becoming warmer and 
closer to observations. The top of Fig_ 3.19b shows that the Control Run 
G-0 has this region almost 5° too cool during the summer. 
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Surface Air Temperature (C) 


Annual Zonally Averaged Temperature 
(Landsurface Only) 


Figure 3.17a 



Annual zonally averaged surface air temperature over landsurface 
areas; observations are from Legates and Willmott (1990) 
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DJF Zonally Averaged Temperature 
(Landsurface Only) 



Latitude 


Figure 3.17b Zonally averaged surface air temperature over landsurface areas 
over landsurface areas for December-January-February; 
observations are from Legates and Willmott (1990) 
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Surface Air Temperature (C) 


JJA Zonally Averaged Temperature 
(Landsurface Only) 



Figure 3.17c 


Zonally averaged surface air temperature over landsurface areas 
over landsurface areas for June— July— August; observations are 
from Legates and Willmott (1990) 


134 



OBSERVED SURFACE AIR TEMPERATURE. JANUARY 






CONTROL (G-0) JAN SURFACE AIR TEMPERATURE (K) 



G — 3 (SP/SO/ST) JAN SURFACE AIR TEMPERATURE (K) 



Figure 3.19a Winter surface air temperatures for Control Run G-0 and 
Run G-3 (space/soil/storm) 
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G-3 (SP/SO/ST) JUL SURFACE AIR TEMPERATURE (K) 



Figure 3.19b Summer surface air temperatures for Control Run G-0 and 
Run G-3 (space/soil/storm) 
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All in all, landsurface temperatures in the GISS GCM are in good 
agreement with observations. We have seen that the modified hydrology of 
the model in Run G-3 has in some areas improved surface air temperatures 
and in some areas made it worse. In the annual zonal mean, however, the 
agreement with data which is already quite good in the current GISS Model II 
has been improved still further with the inclusion of spatial variability. 

4. Precipitation in the GISS GCM 

It has been a recurring theme throughout this section that errors in 
runoff and evaporation can in large part be traced to errors in the 
precipitation generated by the model. Although it is not the point of this 
research to suggest methods of improvement on the precipitation mechanism, it 
is nevertheless instructive to show in more detail the deviation of the model 
precipitation from observed data. 

We have already seen that model precipitation is severely in error over 
many river basins. As shown in Table 3.4a, the model predicts roughly 
2 meters/year rainfall over the Yellow River, yet the observed value is only 
0.68 meters/year. The generally more dry climates of the Yellow, Nile, 
Kolyma, Yukon, Colorado, and Amur river basins are all overpredicted by at 
least 100%. Further, more than half of the 33 selected basins show model 
precipitation to be in error by 50% or more on an annual basis. A more 
complete picture of precipitation can be seen by observing the seasonal 
distributions. 

Figs. 3.20 and 3.21a-b show winter/summer precipitation observations 
(Schutz and Gates, 1971) and those simulated by GCM. The regional high 
rainfall accounts for the overestimated river basin precipitation. Most notable 
is the extreme maxima over Southeast Asia. Generally the 
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Figure 3.20 Global distribution of precipitation for December-February (top) 
and June-August (bottom). 

Observations are from Schutz and Gates (1971) 
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CONTROL (G — 0) DJF PRECIPITATION (MM/DAY) 



G — 3 (SP /30/ST) DJF PRECIPITATION (MM/DAY) 



Figure 3.21a Winter precipitation for Control Run G-0 and Run G 3 
(space/soil/storm) 
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CONTROL (G — 0) JJA PRECIPITATION (MM/DAY) 



G — 3 (SP/SO/ST) JJA PRECIPITATION (MM/DAY) 



Figure 3.21b Summer precipitation for Control Run G-0 and Run G-3 
(space/soil/storm) 
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Figure 3.22 


Gridded Precipitation Comparison: 
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Gridded Precipitation Comparison: 
Tropics (20S - 20N Lat) 



Model Grid Precipitation (mm/day) 
Run G-3 (sp/so/st) 


Gridded precipitation comparison with observations 
1988) over tropics for Control Run G-0 and 
(space/soil/storm) 
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Figure 3.23 
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spatial extent and intensity of the monsoon circulations in the GISS model are 
stronger than that which actually exists. The grid resolution and ocean model 
may be responsible. (Note: Miller and Russell (1990) show that the 4 * 5 

degree resolution version of the GISS GCM has considerably better 
precipitation climatology). 

A precipitation data set from the National Center for Atmospheric 
Research (NCAR) compiled by Spangler and Jenne (1988) was used to 
compare precipitation on a grid square basis. The data set covers the years 
1963 through 1973 and uses more than 3900 different stations over the Globe. 
The data were gridded to match GISS 8° * 10° grid squares by a simple 
"nearest-neighbor" averaging technique (Theissen polygons). Fig. 3.22 shows 
the Control Run and Run G-3 for tropical landsurface grids (20 S to 20 N) 
and Fig. 3.23 shows the same for the Northern Hemisphere (>20 N). This 
latitude splitting shows results similar to the zonally averaged precipitation in 
Fig. 3.13a in the tropics; that is, the mean of the Control Run G-0 is roughly 
equivalent to the mean of observation. On the other hand, the mean of 
Run G-3, upon close inspection, can be seen to be slightly lower than the 
observations. In the Northern Hemisphere, both simulations are slightly higher 
in their means than observations. Thus two independent data sets reveal the 
same diagnosis of rainfall in a rough breakdown of zonal mean values. 
Further, the plots show the severe error of precipitation generation over most 
grids - generally ±50% or greater from observations. 

Finally, we also consider the frequency of rainfall over three selected grid 
squares. Fig. 3.24 shows a map of three U.S. grid squares for which hourly 
precipitation values were recorded during each of the simulations. We define 
these grids as "North Central", "Southeast", and "Southwest" even though 
they cover only part of the areas generally associated with these terms. 
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Figure 3.24 Locations of "North Central U.S.", "Southeast U.S.", and 
"Southwest U.S." grids along with representative stations 
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A data set of hourly precipitation records (spanning 40 years) for major U.S. 
cities (Earthinfo, 1989) was used to compare rainfall frequency. The 
probability of rain in the model is estimated from the frequency of hours for 
which there is either supersaturation or moist-convective rainfall. Similar 
statistics are computed for the observed precipitation for three regions 
corresponding to GISS GCM grids over North America. 

Figs. 3.25a-c show the comparison of the selected GISS GCM grid 
fraction of time with precipitation and the values observed for several 
measurement stations within the grid region. The differences between these 
statistics for Runs G— 0, G— 1, G— 2 and G— 3 are not tested for statistical 
significance since the simulations are only for a few years each. Over the 
three regions, the GISS GCM generally gives between 15 and 20 percent of the 
hours as hours with precipitation. The observations, based on the average of 
seven stations over the "Southeast U.S." region, generally give about 7 percent 
probability of rain in any hour. A value of closer to 6 percent is found for 
the "North Central U.S." region as depicted in Fig. 3.25b. Over the 
"Southwest U.S." region where much of the rainfall is due to moist-convective 
processes, there is a larger range in observed statistic for the five stations 
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Fraction of Time With Precipitation 


Precipitation Frequency Comparison: 
Southeast U.S. Grid 



Figure 3.25a Comparison of Southeast U.S. grid precipitation frequency with 
station data (based on hourly measurements) 
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Fraction of Time With Precipitation 


Precipitation Frequency Comparison: 
North Central U.S. Grid 



Figure 3.25b Comparison of North Central U.S. grid precipitation frequency 
with station data (based on hourly measurements) 
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Fraction of Time With Precipitation 


Precipitation Frequency Comparison: 
Southwest U.S. Grid 



Figure 3.25c Comparison of Southwest U.S. grid precipitation frequency with 
station data (based on hourly measurements) 
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listed in Fig. 3.25c. The fraction of time with precipitation for stations 
enclosed in this region range from 1 percent up to nearly 5 percent. A closer 
look at this region reveals that, due to the spatially-concentrated form in 
which moist-convective rainfall is delivered, a better measure of the grid-wide 
probability of rain is the frequency of hours with precipitation at MI one of 
the stations enclosed within the grid. For regions where rainfall is more 
spatially uniform and the storm coverage is over areas on the same scale as 
the GCM grid or larger (e.g., stratiform rainfall over the "Southeast U.S." or 
"North Central U.S." grids) then the average of precipitation probabilities for 
many stations will be equal to the probability of precipitation at any of the 
many included stations. Fig. 3.26 shows that in the "Southwest U.S." region 
(where rainstorms are of moist-convective origin and are therefore of limited 
spatial extent when compared with the GCM grid area) the fraction of time 
with precipitation is going to be different when computed as the average of 
the statistic for many stations, as opposed to its estimation by considering the 
frequency of times when there is rain at any one of many stations. In 

Fig. 3.26 up to fifteen stations are included for this region. The probability of 
precipitation is computed when groups of 1, 2, 3, •••, 15 stations are 
considered at a time. The figure shows that the statistic is still growing with 
fifteen stations but it is reaching an asymptote. The value of the statistic at 
that asymptote is the true probability of precipitation over the "Southwest 
U.S." grid and it is this statistic that should be compared with the simulation 
results from the GISS GCM. Such comparison is made in Fig. 3.25c with the 
result of considerably improved comparison between the observed and simulated 
probability of precipitation over the "Southwest U.S. grid. 
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Figure 3.26 Rainfall frequency derived from hourly station data within a 
Southwest U.S. GCM grid square for various station-network 
arrays. 

(Data from Earthlnfo, 1989) 
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The GCM-simulated values of regional precipitation are not designed to 
replicate, nor should they be compared with point— measurements. The higher 
order statistics (such as variance, autocorrelation, etc.) of station precipitation 
are thus not useful in validating GCMs. Usually the mean precipitation 
values for grid regions are used in validation studies. As evident in 
Figs. 3.25a-c and Fig. 3.26, the frequency of time-periods with precipitation is 
an additional (non-parametric) measure of the precipitation process which may 
be used to validate models. Fig. 3.26 shows that considerably more stations 
are needed to define robust statistics. Further study is needed to make use of 
this measure in validating GCMs. 

A final observation of hydrologic significance is that potential evaporation 
in the GISS GCM is generally very high (see Figs. A.6a-d). Although outside 
the scope of this work, we direct the reader to Milly (1991) for an analysis of 
alternate formulations of potential evaporation in GCMs. 
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D. Discussion 


We have implemented the landsurface hydrology parameterization of 
Entekhabi and Eagleson (1989a) along with the soil moisture diffusion 

parameterization based on Abramopoulos et al. (1988) into the GISS GCM. 
The results of three new simulations compared against a Control Run and 
observations have been presented. The following is a discussion and evaluation 
of these results. 

The first point to be made is that the new parameterizations are both, 
to a much greater degree than the current GISS II landsurface, physically- 
based. The general problem encountered when improving the model by 
incorporating more realistic schemes, however, is a requirement for large 
amounts of computing time. The times required for our simulations were as 
follows: The Control Run requires roughly one hour CPU time per month of 

simulation on the IBM mainframe computer. Run G-l increased this time by 
approximately 10%, and Runs G-2 and G-3, with the new soil diffusion 

scheme, required about 25% more CPU time than the Control Run. The 
major improvements in the simulations resulted from the inclusion of spatial 
variability. It must be noted that the spatial variability and soil diffusion 
algorithms may be further optimized and their numerics improved in the 
computer code. Thus there is the prospect of much less additional cost if the 
new algorithms are optimally adapted to the GISS GCM. It has been 
demonstrated that with a fractional wetting parameter for moist-convective 
type rainfall, the inclusion of spatial variability results in remarkable 

improvements of the landsurface hydrology in all of the global, zonal, 

continental, and large river basin domains. 
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Another important finding is the strong sensitivity of the correct model 
hydrologic balance on the runoff generation mechanism. Evaporation is only 
reduced by a reduction in available soil moisture through increased runoff. In 
this way, runoff is a major control on actual evaporation in the GISS GCM. 

Precipitation has been seen to be rather poorly represented in the 8 * 10 
degree resolution of the GISS GCM, as comparisons have been made with 
many independent data sets. Russell and Miller (1990) show the 4*5 degree 
resolution version of the GISS GCM has considerably improved precipitation 
climatology. 

The extent to which the landsurface parameterization can affect rainfall 
was seen most clearly in the globally averaged water balance schematics of 
Fig. 3.12. Roughly a 10% reduction resulted over the landsurface due to the 
increased runoff in Run G— 3. However, this is by no means enough to bring 
precipitation in its distribution over the earth into agreement with 

observations. The global average result is good but hides the fact that 
precipitation is actually too low in the tropics and too high in the northern 
hemisphere. In most areas, errors in precipitation create errors in the water 
balance which are impossible to correct with any landsurface parameterization. 
Both the magnitude and temporal and spatial structure of rainfall generation 
are in need of improvement in the atmospheric branch before the landsurface 
water balance will be able to agree with data. Another factor to consider is 
the interannual variability in the model climate and the representativeness of a 

limited number of years of simulation. 

One key area which is of concern when altering the water balance is the 
effect on the heat balance. The diagnostic most available for comparison is 

the surface air temperature. We have seen that in the zonal mean, surface 
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air temperatures have improved over the tropics and have remained unchanged 

elsewhere in Run G-3. This is a reflection of the improvement in evaporation 

over the tropics, giving a more accurate latent heat flux. 

One last point to consider is the lack of data for the two key spatial 

variability parameters. These are 1) the coefficient of variation (cv ) of soil 

s 

moisture over a grid square, and 2) the fractional wetting of moist-convective 

and large scale supersaturation type rainfall. For all simulations, the cv of 

s 

soil moisture was set equal to 1.0. Because of constraints on computer time, 
no GCM simulation tested this sensitivity in the full interactive global case. 
This needs to be done in the future. 

The fractional wetting parameter however was tested since it obviously 
caused strong sensitivity in off-line analyses. This parameter intuitively ought 
to have a dependence on convection, topography, vicinity of oceans, climate 
type, and wind patterns, among other things. We explored only a globally 
uniform setting of this parameter and found strong sensitivity. Further 
consideration ought to be given here for both types of rainfall. One further 
argument in support of a lower value for use in the model is the higher 
frequency of simulated grid rainfall when compared to station observations. If 
grid rainfall occurs more frequently than station data indicate, one could 
reason that the grid precipitation is perhaps accounting for many events each 
of fractional-grid scale. The fractional wetting k needs to be related to the 
fraction (by mass) of the grid air column experiencing convection. This latter 
fraction is explicitly solved for by Kuo-type moist convection schemes. In the 
next generation of the GISS GCM, the land surface hydrology including 
subgrid scale spatial variability may be coupled with the moist convection 
scheme in this manner. The soil moisture diffusion scheme of Abramopoulos 
et al. (1988) also needs to be implemented along with the parallel soil heat 
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diffusion algorithm. Global data sets on soils, vegetation canopy, surface 
topography, roughness, etc. need to be incorporated as well into the improved 
landsurface hydrology schemes. 



Chapter IV 

Conclusions and Recommendations for the giss gcm 
A. Summary of Research Results 

We have implemented improved landsurface hydrological parameterizations 
into the GISS GCM and have analyzed their effects on simulated global 

climate. The statistical-dynamical parameterization of Entekhabi and Eagleson 
(1989a) is used for its critical advantage over current GCM hydrological 
schemes due to the inclusion of subgrid scale spatial variability. The soil 
moisture diffusion parameterization of Abramopoulos et al. (1988) is used to 
provide more realistic deep soil water storage and moisture diffusion. 

Descriptions of these parameterizations are given in Chapter I along with a 

description of an efficient one-dimensional Screening Model with which 

sensitivities of the parameterizations were evaluated prior to implementation in 
the full three-dimensional GISS GCM. 

The results of sensitivity simulations with a three layer soil column in 
the Screening Model are presented in Chapter II. These simulations 
investigate sensitivity to soil storage capacity and heat capacity (i.e., soil layer 
thicknesses), sensitivity to groundwater percolation from the lowest soil layer, 
and sensitivity to transpiration extraction from lower soil layers. It is 

determined that there is strong sensitivity of mean simulated climate to the 
thickness of the top soil layer. The top soil layer thickness also has a strong 
effect on amplitudes of the diurnal heat cycle. For the conditions tested here, 
it is found that increasing the top soil layer thickness created in general, an 
increasingly drier climate while dampening the diurnal heat cycle. Increasing 
the total soil column thickness, while having only a small effect on mean 
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climate, has a significant effect on the annual range of climatic variables. The 
annual ranges of soil moisture and temperature are dampened as the total soil 
column depth increases. Based on these results, the top two soil layer 
thicknesses are held constant over the globe and the lowest soil layer thickness 
adjusted to preserve field capacities of GISS-II upon implementation of the 
Abramopoulos et al. (1988) soil moisture diffusion scheme into the GISS GCM. 

With regard to groundwater percolation, it is found that for conditions of 
the Screening Model simulated climate the sensitivity is significant mainly for 
light textured soils. Simulations with the GISS GCM use a no-flux lower soil 
boundary condition in the absence of terrain and topography data, and 
adequate conceptualization of percolation when considering large areas. 

The Screening Model simulations testing the effects of transpiration from 
lower soil layers show strong sensitivity to the distribution of roots in the soil 
layers. The lower soil layers tend toward excessive drying even with moderate 
fractions of roots in those layers. Also, mean climate is affected most strongly 
by transpiration in the setting of a deep soil column. Lacking a firm data 
base for root distributions in the soil, the choice for distribution used in the 
GISS GCM is based on the Screening Model results and it is set at 85% of 
the roots in the upper soil layer, 10% in the middle layer, and 5% in the 
lowest layer. This distribution produced the most realistic results in the 

Screening Model. 

The preceding discussion now leads to the main focus of this research, 
which is the effects of new landsurface hydrological parameterization in the 
GISS GCM. The parameterizations of Entekhabi and Eagleson (1989a) and 
Abramopoulos et al. (1988) are both physically-based in contrast with the 
current landsurface hydrological parameterization of the GISS GCM. 
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Entekhabi and Eagleson (1989a) account for spatial variability in the key 
parameters of rainfall and soil moisture and incorporate realistic equations of 
infiltration and exfiltration from the soil. Abramopoulos et al. (1988) give a 
finite difference approximation to the governing diffusion equation for soil 
moisture. The additional computation requirements increase the simulation 
run-time by 25% with inclusion of both parameterizations (10% for the spatial 
variability parameterization alone). It is important, however, to note that 
these algorithms may be further optimized in their coding to reduce their cost. 

The results shown in Chapter III.C demonstrate the major improvements 
in the hydrologic balance resulting from the inclusion of spatial variability. 

Results are compared over a wide range of spatial domains (global, zonal, 
continental, and river basins) using a number of data sets and improvements 
are verified on all fronts. Because of the nonlinear response of runoff to soil 

moisture and precipitation intensity when spatial variability is included, the 

strongest changes in hydrologic budgets occur over the tropics. Improvement 
in the hydrologic balance further results in improved heat balance verified 

most distinctly in comparisons of zonal surface air temperature over landsurface 
areas. The soil moisture diffusion scheme of Abramopoulos et al. (1988) is 
necessary for maintaining realistic annual cycles of heat and moisture. 
However, it has smaller effects on the total global hydrologic budget. 

B. The Need for Spatial Variability 

A major finding of this research is the remarkable improvement in the 
landsurface hydrologic balance particularly over the tropics obtained by 
inclusion of spatial variability (see Figures 3.13a-c). The poor agreement of 
the current 8° x 10° GISS-II hydrologic balance can be seen to stem mainly 
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from a lack of runoff generation, giving simulated values much lower than 
those observed over landsurface areas. This allows evaporative fluxes to 
exceed by a large margin the evaporation values derived from observations. 
Based on this research, it seems most likely that the lack of runoff generation 
in the current GISS-II in the tropics is due to the fact that precipitation is 
currently modeled as being uniform over the entire grid square. This results 
in low average intensities which generate far less runoff (if physically-based 
infiltration equations are employed) than would be obtained by spatially 
heterogeneous rainfall, having some areas of concentrated rainfall producing 
larger amounts of runoff and some areas of lesser rainfall intensity producing 
less runoff. This coupled with spatial variability in soil moisture in the 
formulation of Entekhabi and Eagleson (1989a) has been shown here to 
produce results in better agreement with observations. Regardless of the levels 
of detail which may be pursued in modeling landsurface hydrology, without the 
element of spatial variability it seems unlikely that the global hydrologic 
balance will be represented adequately in GCMs. 

C. Future Research 

In order to further improve the landsurface hydrological parameterization 
of GCMs several avenues of research need to be explored. First, it is 
determined that there is very strong sensitivity in the spatial variability 
parameterization to the rainfall fractional wetting parameter k. As mentioned, 
this parameter ought to have a dependence on the air column convection, 
topography, seasonality, and prevailing climate among other factors. While 
simulations here simply used one value of k for moist-convective rainfall and 
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one for supersaturation rainfall, alternate formulations need to be investigated. 
The same is true for the coefficient of variation of soil moisture cv s . 

Groundwater percolation is treated only in a very simple fashion in this 
work. In nature, groundwater percolation can be a major component in the 
hydrologic cycle and as such it needs to be investigated and its influence 
quantified in the context of GCMs. 

In the sensitivity experiments with the One-Dimensional Screening Model, 
it is found that the top soil layer and the total soil depth determine the mean 
climate and the amplitude of its diurnal and seasonal cycles. There is a 
critical need to clarify this sensitivity further and develop objective methods 
by which the soil column may be discretized. There is also the need to 
search sources of data for defining this important lower boundary. 

The simulations performed in the GISS GCM here have been only of a 
relatively short duration (maximum of five years), using fixed sea surface 
temperatures. Longer duration simulations may reveal trends and statistical 
measures by which the model could be further analyzed. An interactive ocean 
component ought to be used as well. 

In Chapter III, a new measure is defined to be used in validating the 
precipitation process of the GCM model climate. Since the second-order 
properties (variance, covariance) of GCM— produced rainfall cannot be compared 
with that resulting from weather at an observation point, the new statistic 
defined here will be especially useful since it measures the structure of 
variability of the modeled rainfall without relying on second-order statistics. 
Together with the mean, the probability of regional precipitation may be 
estimated from observed data and employed in validating the precipitation 
climatology of GCM model climates. 
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Finally, because of their strong impact on the landsurface hydrologic 
budget, the GCM generated rainfall distributions and potential evaporation 
mechanism require a more thorough examination. Without the proper 
potential evaporation forcing and rainfall generation, the landsurface will not 
be able to partition these forcings accurately in order to represent the 
hydrologic balance. 
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Appendix 


A. Seasonal Fields of Simulated Hydrologic Diagnostics 

The following figures (A.l.a through A.6.d) show the seasonal fields of 
simulated precipitation, evaporation, surface runoff, relative surface soil 
saturation, surface air temperature, and potential evaporation for the 
simulations G-0 through G-3 (simulations with the GISS GCM). The 
presentation format gives the Control Run (G— 0) first in actual diagnostic 
values followed by fields of the differences between each of the simulations 
G-l through G-3 and the Control Run. 
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B. List of Appendix Figures 


Figure A.l.a 
Figure A.l.b 

Figure A.l.c 

Figure A.l.d 

Figure A.2.a 
Figure A.2.b 

Figure A.2.c 

Figure A.2.d 

Figure A.3.a 
Figure A.3.b 

Figure A.3.C 

Figure A.3.d 


Precipitation Fields 

Seasonal precipitation fields (mm/day) of Control Run (G-0) 

Seasonal precipitation fields (mm/day) of G-l (presented as 
[(G— 1) - Control]) 

Seasonal precipitation fields (mm/day) of Run G-2 (presented 
as [(G-2) - Control]) 

Seasonal precipitation fields (mm/day) of Run G-3 (presented 
as [(G-3) - Control]) 


Evaporation Fields 

Seasonal evaporation fields (mm/day) of Control Run (G-0) 

Seasonal evaporation fields (mm/day) of Run G-l (presented 
as [(G-l) - Control]) 

Seasonal evaporation fields (mm/day) of Run G-2 (presented 
as [(G-2) - Control]) 

Seasonal evaporation fields (mm/day) of Run G-3 (presented 
as [(G-3) - Control]) 


Runoff Fields 

Seasonal runoff fields (mm/day) of Control Run (G-0) 

Seasonal runoff fields (mm/day) of Run G-l (presented as 

[(G-l) - Control]) 

Seasonal runoff fields (mm/day) of Run G-2 (presented as 

[(G-2) - Control]) 

Seasonal runoff fields (mm/day) of Run G-3 (presented as 

[(G-3) - Control]) 
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Figure A.4.a 
Figure A.4.b 
Figure A.4.c 
Figure A.4.d 

Figure A.5.a 
Figure A.5.b 
Figure A.5.C 
Figure A.5.d 

Figure A.6.a 
Figure A.6.b 
Figure A.6.C 
Figure A.6.d 


Soil Saturation Fields 

Seasonal relative surface soil moisture fields (%) of Control 
Run (G— 0) 

Seasonal relative surface soil moisture fields (%) of Run G-l 
(presented as [(G-l)- Control]) 

Seasonal relative surface soil moisture fields (%) of Run G-2 
(presented as [(G-2)- Control]) 

Seasonal relative surface soil moisture fields (%) of Run G— 3 
(presented as [(G-3)- Control]) 


Surface Air Temperature Fields 

Seasonal surface air temperature fields (C) of Control Run 
(G-0) 

Seasonal surface air temperature fields (C) of Run G-l 

(presented as [(G-l) - Control]) 

Seasonal surface air temperature fields (C) of Run G-2 

(presented as [(G-2) - Control]) 

Seasonal surface air temperature fields (C) of Run G-3 

(presented as [(G-3) - Control]) 


Potential Evaporation Fields 

Seasonal potential evaporation fields (mm/day) of Control Run 

(G-0) 

Seasonal potential evaporation fields (mm/day) of Run G-l 

(presented as [(G-l) - Control]) 

Seasonal potential evaporation fields (mm/day) of Run G-2 

(presented as [(G-2) - Control]) 

Seasonal potential evaporation fields (mm/day) of Run G-3 

(presented as [(G-3) - Control]) 
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Figure A.l.a Seasonal precipitation fields (mm/day) of Control Run (G-0) 
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[(G — 1) — CONTROL] DJF PRECIPITATION (MM/DAY) 


[(G — 2) — CONTROL] DJF PRECIPITATION (MM/DAY) 


[(G — 3) -CONTROL] DJF PRECIPITATION (MM/DAY) 



[(G— 3) -CONTROL] JJA PRECIPITATION (MM/DAY) 



Figure A.l.d Seasonal precipitation fields (mm/day) of Run G-3 (presented 
as [(G-3) - Control]) 
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[(G — 1) -CONTROL] DJF EVAPORATION (MM/DAY) 



[( G — 1.) — CONTROL] JJA EVAPORATION (MM/DAY) 



Pieure A.2.b Seasonal evaporation fields (mm/day) of Run G-l (presented 
as [(G-l) - Control]) 
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Figure A.2.c Seasonal evaporation fields (mm/day) of Run G-2 (presented 
as [(G-2) - Control]) 



[(G- 3) -CONTROL] DJF EVAPORATION (MM/DAY) 



[(G-3) -CONTROL] JJA EVAPORATION (MM/DAY) 



CONTROL DJF SURFACE RUNOFF (MM/DAY) 



CONTROL JJA SURFACE RUNOFF (MM/DAY) 



[(G-l)-CONTROL] DJF RUNOFF (MM/DAY) 


[(G- 2) — CONTROL] DJF RUNOFF (MM/DAY) 


[(G- 3) — CONTROL] DJF RUNOFF (MM/DAY) 


CONTROL DJF RELATIVE SURF SOIL MOISTURE (%) 



CONTROL JJA RELATIVE SURF SOIL MOISTURE (%) 





[(G— 1)— CONTROL] DJF SOIL MOISTURE (%) 



[(G— 1) - CONTROL] JJA SOIL MOISTURE (%) 



Firare A.4.b Seasonal relative surface soil moisture fields (%) of Run G 1 
(presented as [(G-l)- Control]) 
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Figure A.4.c Seasonal relative surface soil moisture fields (%) of Run G-2 
(presented as [(G-2)- Control]) 
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[(G-l)-CONTROL] DJF SURFACE AIR TEMP (C) 


[(G— 2) -CONTROL] DJF SURFACE AIR TEMP (C) 



[(G — 2)— CONTROL] JJA SURFACE AIR TEMP (C) 



Figure A.5.C Seasonal surface air temperature fields (C) of Run G-2 
(presented as [(G-2) - Control]) 
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[(G— 3) -CONTROL] DJF SURFACE AIR TEMP (C) 


CONTROL DJF POTENTIAL EVAP (MM/DAY) 


[(G— 1) — CONTROL] DJF POTENTIAL EVAP (MM/DAY) 



[(G-l)-CONTROL] JJA POTENTIAL EVAP (MM/DAY) 



Figure A.6.b Seasonal potential evaporation fields (mm/day) of Run G 1 
(presented as [(G-l) - Control]) 
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[(G-2) -CONTROL] DJF POTENTIAL EVAP (MM/DAY) 



[(G — 2) -CONTROL] JJA POTENTIAL EVAP (MM/DAY) 



Figure A.6.C Seasonal potential evaporation fields (mm/day) of Run G-2 
(presented as [(G-2) - Control]) 
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[(G — 3) — CONTROL] DJF POTENTIAL EVAP (MM/DAY) 



